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This is a Semi-Annual Siams Report on research conducted between 22 September 1990 
an.,21 Match 1991 under NASA Grant NAG 5-814, entitled "The Interpretation of Crustal 
Dynamics Data in Terms of Plate Motions and Reghmal Defonnadon near Plate Boundane, 

This gran, supports the research of one Investigator (S. C. Solomon), one Research Staff (R. 
Reilinger), and two Ph. D. students (A. F. Sheehan andC. I. Wolfe) on behalf of the NASA 

Geodynamics and Crustal Dynamics Programs. 

The focus of the research has been in mo brcal arcas during the most recent fr month 

period- (1) the nature arc! dynamics of time-dependent deformation and stress along major 
seismic zones, and (2) the namrc of long-waveiengti, oceanic geoid anomalies in terms of 
lateral variations in upper mantle temperature and composition. The principal ftndrngs of our 
research to date are described in the accompanying appends, m firs, two and the fourth are 
preprints of papers recently submined for publication, and the tithd is the absnac. of a recently 

completed Ph.D. thesis supported by this project. 
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GPS Measurements of Strain Accumulation 
Across the Imperial Valley, California: 
1986-1080 


Abstract 


GPS data collected in southern California from 1086 to 1989 indicate 
considerable strain accumulation across the Imperial Valley. Displacements 
are computed at 29 stations in and near the valley from 1986 to 1988, and at 
11 sites from 1988 to 1989. The earlier measurements indicate 6.9 ± 1.0 
cm/yr right-lateral differential velocity across the valley, although the data 
are heavily influenced by the 1987 Superstition Hills earthquake sequence. 
Some measurements, especially the east-west trending displacements, are 
suspect for large errors. The 1988-1989 GPS displacements are best modeled 
by 5.2 ± 0.9 cm/yr of plate^boundary deformation, but rates calculated from 
conventional geodetic measurements (3.4-4.3 cm/yr) fit the data nearly as 
well. There is evidence from GPS and VEBI observations that the present slip 
rate along the southern San Andreas fault is smaller than the long-term 
geologic estimate, suggesting a lower earthquake potential than is currently 
assumed. The Imperial Valley GPS sites form part of a larger network of 183 


stations spanning an entire cross-section of southern California and northern 
Mexico. Once data from a recent 1990 campaign are fully analyzed and 
integrated with the previous measurements, the strain distribution across the 
San Andreas, San Jacinto, and Elsinore faults will be well established. 
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5.1 Introduction 

The Global Positioning System (GPS) is rapidly becoming one of the most 
important tools to study tectonic deformation. Signals from earth-orbiting 
NAVSTAR satellites (NAVigation Satellite Time And Ranging) are inverted 
to obtain 3-dimensional coordinates of geodetic monuments with high 
precision. For crustal deformation studies, the relative position (or baseline) 
between stations is often measured. Under optimal conditions, the typical 
accuracy for a 50 km baseline is about 1 cm in the horizontal and 3 cm in the 
vertical [e.g., Dams et ai, 1989). The accuracy is significantly degraded from 
poor observing conditions. GPS measurements are used to monitor secular 
deformation such as that associated with plate motion, or to record rapid 
strain fluctuations such as those due to seismic and volcanic activity. 

A prime location for GPS studies is the Imperial Valley of southern-most 
California (Figure 5.1). The valley is one of the most tectonically active 
regions in the state and has been the site of several large earthquakes. In 
fact, GPS monitoring was initiated in 1986 with resurveys in 1988 and 1989. 
GPS station displacements from 1986 to 1988 have been discussed by Larsen 
et ai (1990). These measurements illustrate the effect of the 1987 Superstition 
Hills earthquake sequence, as well as nonseismic displacement components 
attributed to interplate deformation. In the present study we incorporate the 
1989 GPS observations. Station displacements from 1986-1988 and 1988-1989 
are computed and these movements are related to the relative motion between 

the North American and Pacific plates. 


5.2 Seismicity and Tectonics 
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Figure 5.1: Major faults and seismicity from 1832 to 1980 (Caltech 

Catalog) in the Imperial Valley. Large earthquakes are shown as stars. 

The Brawley Seismic Zone is the region of anomalously high activity 

between the Imperial and San Andreas faults. Major earthquakes include 

the 1940 and 1979 events along the Imperial fault, the 1954 and 1968 

events along the San Jacinto fault, and the 1987 Superstition Hills 

earthquake sequence along the Superstition Hills and Elmore Ranch 
faults. 
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The Imperial Valley (Figure 5.1) is a complex transition zone between 
crustal spreading in the Gulf of California and right-lateral transform motion 
along the San Andreas fault [Lommtz cl cl, 1870; Eldcre cl of., 1072|. The 
valley is 4-5 million years old and has been filled by up to 5 km of late 
Cenozoic sediments [Ari. cl cl., 1082). The structural axis of the valley and 
Us major fault systems trend to the northwest, roughly parallel to the 
Pacific-North American plate motion. A significant fraction of the relative 
plate displacement may be accommodated across the valley. 

The Imperial Valley is one of the most seismically active regions of 
California with much of the activity occurring along the Imperial fault and in 
the Brawley Seismic Zone [Johnson and Hill, 1982). Several large earthquakes 
have occurred in and near the Imperial Valley since 1040. The Imperial fault 
ruptured with a M s 7.1 event in 1040 and a M L 6.6 event in 1970 f IJ s G S 
1982). Segments of the San Jacinto fault system broke with a M L 6.2 
earthquake in 1954 (Clark segment) and the 1968 M L 6.5 Borrego Mountain 
event (Coyote Creek segment). The most GPS relevant episode of seismic 
activity occurred recently along the Superstition Hills segment of the San 
Jacinto fault system |e.g„ Me, {circle cl cl., 1080). On November 24, 1087, a 
large (Ms 6.2) earthquake occurred along a northeast trending seismic 

lineament, and was followed 12 hours by a larger event (M s 6.6) along the 
Superstition Hills fault. 

Conventional geodetic measurements indicate significant displacements 
across the Imperial Valley, inferred to represent interplate deformation. 
Triangulation data averaged from 1941 to 1986 suggest 4.3 cm/yr plate- 
boundary movement oriented N40’ W [Snsy end Drew, 1988). The observed 
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deformation is time dependent, with rates of 6.1, 2.1, and 4.5 cm/yr for the 
intervals 1941-1954, 1954-1967, and 1967-1979, respectively. The high 
velocity for the earliest period supports the hypothesis of northwestward 
strain migration following the 1940 earthquake {Thatcher, 1979; Re ilmger, 
1084]- Furthermore, the computed station displacements indicate that north 
of the Imperial fault interplate deformation is distributed over a tone at least 
50 km wide whereas to the south interplate deformation is concentrated 
within a 20 km wide band centered along the Imperial fault. Trilateration 
measurements made by the U.S. Geological Survey from 1972 to 1987 
(Frescuit el ai, 1987a; Prescott et ai, 1987b| indicate 3.45 cm/yr differential 
displacement oriented roughly N40'W between stations on opposite sides of 
the valley. Unlike the increased rate following the 1940 earthquake, these 
measurements reveal no significant change in station velocity after the 1979 
event \Savagz ct a/., 1986], 

New global plate models (NUVEL-1) [DcMcts ct ai, 1987; DcMcts ct ai, 
1990] predict the Pacific-North American relative velocity at Imperial Valley 
coordinates (longitude 115.5‘W, latitude 33.0 *N) averaged over the last 
several million years to be 4.7 cm/yr oriented N39.6*W. VLBI observations 
during the 1980’s suggest a similar present-day rate [e.g., Clark et ai, 1987; 
Kroger et ai, 1987). A significant fraction of this motion may be distributed 
along faults in and near the Imperial Valley. 

5.3 GPS Observations 

The data presented here were obtained in a series of GPS field campaigns 
in 1986, 1988, and 1989. Here we provide summaries regarding the collection 



-270- 


and processing of the 1986 and 1988 surveys, discussed in more detail by 
Larsen et al. [I990j, and present a more detailed account of the most recent 
campaign. In all, a total of 32 Imperial Valley stations have been occupied 
more than once between 1986 and 1989 (Table 5.1). 

The National Geodetic Survey (NGS) began GPS observations in southern 
California with a 54 station network in 1986; 42 stations are located in or 
near the Imperial Valley (Figure 5.2). TI-4100 GPS receivers were used for all 
data collection. Each of the 20 days of observation are processed 
independently utilizing the GPS22 software developed at the NGS, with 
satellite orbit parameters provided by the NSWC (Naval Surface Weapons 
Center). The day to day solutions are integrated into one set of station 
coordinates with the geodetic adjustment program DYNAP (DYNamic 
Adjustment Program) [ Drew and Snay , 1989). Due to insufficient coverage 
and poor data quality during the 1986 survey, uncertainties are suggested to 
be approximately 1 ppm (parts per million) [Neugcbauer, 1988). 

During late February and early March, 1988, university GPS crews 
(UNAVCO) occupied 19 sites in the Imperial Valley, including 15 marks 
observed in 1986 (Figure 5.2). The following month, the NGS returned to the 
Imperial Valley and reoccupied 21 previously established monuments. TI- 
4100 receivers were used for all measurements. Data from both surveys are 
processed independently with the Bernese GPS software package from the 
University of Bern in Switzerland. For each campaign, the data are combined 
into one multiday solution. Satellite orbits are improved with the aid of 
fiducial observations from Mojave (California), Westford (Massachusetts), and 
Richmond (Florida), made as part of the Cooperative International GPS 
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Table 5.1 Reoccupied Stations (1986-1989) 


Black Butte NCMN 1982 BLAC 

Brawley 2 rm 5 BRAW 

Calexico 1954 CALE 

Calipatria 2 CALI 

Coach COAC 

College 1967 COLL 

El Centro 2 1959 ELEC 

Frink 1934 ERIN 

GLO Comew 1934 GLOC 

Hamar 2 1967 HAMA 

Holt 1924 HOLT 

Holtville (Alt) 1934 HLTV 

Imp 1934 IMP! 

Imperial 1934 IMPE 

Junction JUNC 

Kane 1939 KANE 

L 589 1967 L589 

Mack 2 1967 bm reset MACK 

Mello 3 1967 MELL 

Monument Peak NCMN 1983 MONU 

OcotiUo NCMN 1982 OCOT 

Ocotillo 1935 OCTI 

Offset 217 0217 

Offset 224 0224 

Offset 227 0227 

Orient 1939 ORIE 

Pinyon Flat PINY 

Sandy Beach SANl 

T 1226 T122 

Tamarisk 3 1967 TAMA 


• 

-115.6111 

33.1964 

-72.91 

• 

-115.7198 

33.6638 

490.01 


-115.5434 

32.9773 

-66.78 


-115.5064 

32.6645 

-34.07 

• 

-115.5088 

33.1690 

-89.03 


-115.4070 

33.1962 

-8.40 

• 

-115.5024 

32.8269 

-54.63 


-115.5622 

32.7846 

-45.76 

• 

-115.6470 

33.3603 

-85.39 


-115.2465 

32.8396 

-15.46 


-115.5007 

33.0375 

-79.80 


-115.3963 

32.7814 

-36.79 


-115.3821 

32.8084 

-38.77 


-115.5698 

32.8982 

-59.27 


-115.5788 

32.8439 

-51.85 

• 

-115.0619 

32.7092 

8.35 

• 

-115.8237 

33.0614 

10.05 

• 

-115.7611 

32.9506 

13.62 


-115.1441 

32.7288 

1.20 


-115.4653 

32.7961 

-47.00 

• 

-116.4228 

32.8918 

1839.41 

• 

-115.7962 

32.7901 

-36.43 

• 

-116.0017 

32.7338 

111.33 


-115.3131 

32.6782 

-17.16 


-115.7055 

32.8493 

50.00 


-115.8173 

32.6414 

70.64 

• 

-115.4064 

32.9168 

-61.36 


-116.4588 

33.6093 

1235.88 

• 

-115.8344 

33.1929 

-99.92 


-114.8050 

32.8583 

141.19 


-115.4783 

32.8829 

-68.19 


Station Name 


Abbr. 


Occupation Longitude Latitude Elevation 
1986 1988 1989 1 (s) 


Acute 1934 


ACUT 


-115.6093 33.0300 -80.91 
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Figure 6.2: GPS stations surveyed in 1986 and 1988. The 1986 
campaign was conducted by the National Geodetic Survey and included 
42 stations in and near the Imperial Valley. The 1988 observations 
consisted of two campaigns, the first by university groups in 
February/March and the second by the National Geodetic Survey in 
March/April. A total of 32 stations were occupied in 1988, of which 29 

were repeat measurements from 1986. Stations mentioned in text are 
indicated. 



L6 


.5 











- 274 - 


Network (CIGNET) [Chin, 1987). Since orbit improvement techniques are 
used, the (horizontal) precision for each survey is about ~ 0.03 ppm (sub- 
centimeter) [e.g., Davis tt al., 1989; Dong and Bock, 1989]. The Cartesian 
coordinate differences from the university and NGS surveys are adjusted by 
least squares to obtain station positions for 1988. 

During March, 1989, university groups occupied 28 geodetic marks in the 
vicinity of the Imperial Valley, 19 of which were previously surveyed in 1986 
or 1988 (Table 5.2, Figure 5.3). Several new marks were established north of 
the Salton Sea in the Coachella Valley. While most data were collected with 
TI-4100 GPS instruments, this campaign differed from previous surveys in 
that Trimble-4000SD receivers were used at some sites. The field experiment 
was conducted at a time of anomalously high solar flare activity which 
created large ionospheric disturbances [Jackson ct al., 1989). The ionosphere 
creates a frequency dependent delay for the GPS multi-signal structure, 
composed of two carrier phase transmissions at 1575.42 MHz (Ll) and 1227.60 
MHz (L2) [e.g., Kmg ct al., 1985). For dual frequency observations (Ll and 
L2) the ionospheric contribution (error) is removed by an appropriate 
combination of the two phase observables. However, if only single frequency 
measurements are available (either intentionally or due to poor observing 
conditions), the positioning accuracy on all but the shortest baselines will be 
seriously degraded. The 1989 phase observations contain a disproportionate 
number of cycle slips and data gaps, presumably due to the poor ionospheric 
conditions. The TI-4100 instruments generally collected both the Ll and L2 
phase signals, so the ionospheric effect could be eliminated. The Trimble 
4000SD receivers, however, experienced significant difficulty maintaining 
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Table 5.2 Imperial Valley GPS Occupation History (1989) 


) ay Stations 

M OCOT PINY' MONU BLAC ROBO EXTR 1 COAC' 

68 OCOT PINY 1 MONU BLAC ROBO EXTR 1 COAC 1 

67 OCOT PINY 1 MONU BLAC VARN SANl 1 

68 OCOT TRAN 1 FRIN COCH VARN SANl 1 

69 OCOT TRAN 1 FRIN DUNP COCH GLOC 1 VORO 1 

70 OCOT TRAN 1 ALAM DUNP ROBO GLOC 1 VORO 1 

72 OCOT KANE ALAM L589 TAMA* GLOC 1 

73 OCOT KANE 0217 1 L589 TAMA 1 ACUT 1 

74 OCOT BORD 0217 1 ORIE SIPH ACUT 1 

75 OCOT BORD 1 CALI ORIE TAMA* SIPH 

78 OCOT JUNC COLL 1 HAMA 1 TAMA 1 SIPH SANl 

77 JUNC COLL ORIE 1 TAMA 1 GLOC 1 SANl 


Day is Julian day of year 
1 Trimble 4000SD 
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Figure 5.3: Imperial Valley GPS stations surveyed in 1989. TI-4100 
GPS receivers (triangles) were used at most sites. Trimble 4000sd 
receivers (open circles) were also used. Thirty sites were occupied; 10 for 
the first time. Due to very poor ionospheric conditions, data collected 
with the Trimble 4000SD receivers are not discussed here. 



32.5 
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phase-lock on the L2 frequency (it is found that newer Trimble models, 
specifically the 4000SDT, are not as susceptible to solar activity). In fact, 
between 30-60 percent of the L2 data (Trimble 4000SD) was lost. It is 
unlikely that the centimeter level accuracy required for this study could be 
achieved solely with the Ll frequency. Therefore, data collected with Trimble 
4000SD instruments are not considered, although we are currently working on 
schemes to utilize these measurements through ionospheric modeling 
constrained by the dual frequency TI-4100 data. Continental fiducial phase 
observations from the CIGNET tracking sites were either nonexistent or of 
extremely poor quality, presumably due to the unfortunate ionospheric 
conditions. We were therefore not able to apply orbit improvement 
techniques so a multiday solution is obtained with the Bernese software 
utilizing the broadcast orbits. Positioning errors from the broadcast 
ephemerides are believed to be 0.1-1. 0 ppm. 


5*4 GPS Displacements 

GPS vector displacements for the intervals 1986-1988 and 1988-1989 are 
shown in Figures 5.4 and 5.5, respectively. AH measurements are made 
relative to station OCTI. Formal estimates of GPS uncertainty almost 
always underestimate variances derived from repeatability studies. We 
attempt to define more realistic errors by multiplying the structure of the 
formal covariance matrix calculated with the GPS solution by an estimated 
variance factor, which scales as the average baseline length. For the 1986- 
1988 displacements, we assume a variance factor so that the average baseline 
error is 5 cm, or 1 ppm for a 50 km baseline. The large uncertainty is due to 
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the poor quality of the 1986 data. For the 1988-1989 displacements, the 
average baseline error is assumed 2 cm, or 0.4 ppm for a 50 km baseline. The 
largest component of uncertainty is attributed to the broadcast orbits used 
for the 1989 solution. Larsen et al. [1990] found 1-3 cm discrepancies between 
broadcast and improved-orbit solutions in a similar sized network spanning 
the Santa Barbara channel. This approach, albeit somewhat ad hoc, allows 
for self consistent relative errors and it illustrates the much larger 
uncertainties in the east-west direction (~ 4 times larger than the north-south 
uncertainties). This distortion is primarily due to the north-south ground 
track of the satellite orbits, which significantly improves solution constraint 
along this orientation. 

Displacements for the 1986-1988 interval (Figure 5.4) are complicated by 
the 1987 Superstition Hills earthquake sequence, as well as large measurement 
uncertainties. Effects of the earthquake sequence are clearly demonstrated in 
the GPS vectors; displacements at KANE and L589 approach 0.5 meters. 
Estimates of fault rupture suggest 10 stations near the seismic rupture zone 
moved at least 5 cm [Larsen et al., 1990]. The displacements are consistent 
with right-lateral slip along the Superstition Hills fault and left-lateral slip 
along the Elmore Ranch fault. Still, there is a considerable component of 
southeast trending movement which can not be explained as seismic 
deformation or measurement uncertainty. The relative displacement between 
stations on opposite sides of the valley averages 5-6 cm/yr. We take this 
motion to represent continuous strain accumulation due to plate motion. 

The 1988-1989 station displacements clearly demonstrate the right-lateral 
southeast trending movement across the GPS network. Stations furthest to 
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Figure 5.4: GPS station displacements for the interval 1986-1088 (1.8 
years). All measurements are made relative to station OCTI. Errors are 
determined by multiplying the formal uncertainties from the GPS 
solution by a variance factor so that the average baseline error scales as 1 
ppm. The east-west uncertainties are about 4 times larger then the 
north-south. Seismically induced displacements from the 1987 
Superstition Hills earthquake sequence are most apparent at stations 
KANE and L589. The large non-seismic displacements are assumed to 
represent relative motion between the Pacific and North American plates, 
which is concentrated across the valley. 
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Figure 5.5: GPS station displacements for the interval 1988-1989 (1.0 
years). All measurements are made relative to station OCTI. Errors are 
determined by multiplying the formal uncertainties from the GPS 
solution by a variance factor so that the average baseline error scales as 
0.5 ppm. Stations to the northeast moved about 5 cm southwest relative 
to stations on the other side of the valley. 
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the northeast are displaced approximately 5 cm to the southeast relative to 
sites on the other side of the valley, although some of the observed motion 
(e.g„ BLAC) may be distortion from the larger east-west uncertainties. In 
addition. Figure 5.5 demonstrates how easily GPS can monitor tectonic 
deformation even over time scales as short as 1 year. 


The 1986-1888 and 1988-1989 displacements are decomposed into their 
north-south and east-west components, corresponding to the directions of 
minimum and maximum error, respectively (Figure 5.6 and 5.7). Each 


component is plotted as the distance from OCTI on a cross section trending 

NSO-E. perpendicular to the predicted plate motion orientation (N40‘W) 

(Figure 5.8). Simple dislocation theory |e.g„ AfonsmAc end Smylic, 1971) is 

used to remove the effect of the 1987 Superstition Hills earthquake sequence 

from the observed 1986-1988 displacement field, following fault models 

suggested by iersen ct at. (1990) (approximately 109 cm right-lateral slip on 

the Superstition Hills fault and 45 cm left-lateral slip on the Elmore Ranch 
fault). 


Decomposing the vector displacements into geographic components tends 
to separate the uncertainties which are magnified along the longitudinal 
direction. The north-south movements clearly exhibit right-lateral 
displacement for both intervals; stations to the northeast display southerly 
offsets relative to sites on the other side of the valley. The magnitude of the 
displacement is roughly proportional to the time interval spanned by the 
measurements, suggesting continuous strain accumulation. Stations which 
display the largest scatter in Figure 5.6 are for the most part those sites where 
the applied seismic correction is greater than 4 cm (open circles). This may 
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indicate additional fault complexity not accounted for by the dislocation 
model used to remove the effects of the 1987 earthquake. 

The east-west movements for the 1986-1988 interval exhibit large scatter 
with no discernible trend across the valley. This is invariant of the seismic 
correction, so the scatter is not explained simply as unmodeled effects from 
the Superstition Hills earthquakes. Presumably, the large deviations are due 
to fairly significant east-west oriented errors from the 1986 survey. This may 
explain the anomalous vector displacements observed in Figure 5.4, especially 
noticeable for those sites near the border east of the Imperial fault. On the 
other hand, the 1988-1989 displacements clearly show large east directed 
movements for stations furthest to the northeast, consistent with southeast 
trending deformation across the valley. 

5.5 Discussion 
Deformation across valley 

The 1986-1988 measurements are concentrated along the Imperial fault 
(Figure 5.4). The nonseismic displacements reveal a sharp boundary 15-20 
km wide between deformation on either side of the valley (Figure 5.6). This 
suggests that in the southern half of the Imperial Valley strain is being 
accommodated exclusively along the Imperial fault. The 1988-1989 
measurements are distributed more uniformly throughout the region (Figure 
5.5), and indicate a broader strain-transition zone (Figure 5.7). This implies 
that deformation may be occurring along several structures to the north, 
including the San Andreas, San Jacinto, and Elsinore faults. The same 
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Figure 5.8: The north-south and east-west displacement components for 
the 1986-1988 interval. All distances are relative to OCTI on a cross 
section trending N50*E, perpendicular to the plate motion (see Figure 
5.8). The effect of the 1987 Superstition Hills earthquake sequence is 
removed. Open circles indicate stations where the seismic correction is 
greater than 4 cm. The north-south offset between stations on opposite 
sides of the valley is 8.1 cm. The large scatter for the east-west 
components are presumably due to errors in the 1986 survey. The 
average uncertainty for each displacement component is shown. 
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Figure 5.7: The north-south and east-west displacement components for 
the 1988-1989 interval. All distances are relative to OCTI on a cross 
section trending N50 * E, perpendicular to the plate motion direction. The 
data are best-fit by 5.2 cm/yr displacement across the valley (solid line), 
although a rate of 3.4 cm/yr fit the data nearly as well (dashed line). 
The average uncertainty for each displacement component is shown. 
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Figure 6.8: Shear plane (10 km depth) used to model the 1988-1989 
displacements (shaded band); cross section used in Figures 5.6 and 5.7; 
and stations surveyed at least twice between 1986 and 1989. 
Considerable strain is observed across the GPS network, which is 

attributed to plate-boundary deformation between the North American 
and Pacific plates. 
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pattern is observed in the conventional geodetic measurements, which indicate 
concentrated strain in a narrow 20 km wide zone about the Imperial fault, 
and diffuse deformation of at least 50 km wide to the north [Snay and Drew, 
1988; Prescott et ai, 1987b]. 

The 1986-1988 displacements (Figure 5.4) have been modeled by Larsen et 
al. [1990]. The station movements are consistent with right-lateral slip along 
the Superstition Hills fault and left-lateral offset along the Elmore ranch fault. 
The nonseismic residuals indicate remaining deformation across the valley, 
evidenced by the 8.1 ± 1.3 cm offset in the north-south displacement 
component (Figure 5.6). The differential movement is calculated by linearly 
fitting those data furthest to the southwest and northeast. The data errors 
are increased by 0.33 times the estimated seismic displacements, giving less 
weight to those stations most affected by the 1987 earthquakes. The east- 
west components are not used due of the large data scatter. The observed 
movement is taken to represent plate-boundary deformation due to the 
relative motion between the North American and Pacific plates. The 
calculated north-south differential displacement averaged over the 1.8 year 
observation interval, is equivalent to 5.9 ± 1.0 cm/yr right-lateral movement 
oriented N40 * W, assuming a uniform velocity field parallel to the direction of 
plate motion. Conventional geodetic data are consistent with the assumed 
orientation [ Prescott et ai, 1987a; Snay et al., 1988]. 

The differential movement across the Imperial Valley is clearly visible in 
the 1988-1989 displacements (Figure 5.5 and 5.7), and is observed in both the 
north-south and east-west components. The movement is smaller than the 
earlier interval because of the shorter observation period (1.0 years). 
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However, the most recent measurements are not influenced by seismic activity 
and contain smaller experimental error. Since the 1988-1989 station 
displacements are more uniformly distributed across the valley, it is difficult 
to constrain an absolute differential ofiset. Instead, the measurements are 
modeled assuming a semi-infinite right-lateral shear plane at depth 
representing the Pacific- North American plate margin (Figure 5.8). The plane 
is oriented N40*W about coordinates 32.796 *N, 115.454 # W, almost 

congruent with the Imperial fault and the axis of the valley. The upper depth 
is constrained at 10 km and uniform slip is assumed over the entire shear 
boundary. Snay and Drew [1988] incorporate a similar model to explain 
triangulation observations from 1941 to 1986, but allow additional slip along 
the Imperial fault necessitated by the detailed station coverage in this region. 
More complex models assuming distributed offset along the Imperial, San 
Andreas, San Jacinto, and Elsinore faults, and within the Brawley Seismic 
Zone have been used to explain other geodetic measurements in the valley 
[e.g., Savage et ai , 1979]. The measurements presented here are not of 
sufficient resolution or accuracy (due to the short time coverage) to warrant 
such detail. The 1988-1989 GPS displacement vectors are best constrained by 
5.2 ± 0.9 cm/yr plate-boundary deformation. The best-fit solution to the 
observed GPS movements is shown in Figure 5.7. Additional solutions are 
obtained by varying the depth to the upper boundary of the shear plane from 
5 to 15 km. The calculated displacement rates range from 4.4 (5 km) to 6.0 
cm/yr (15 km), while the minimum residual solution is obtained at 10 km 
depth (5.2 cm/yr). 

Because the 1988-1989 displacements are not affected by seismic 
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deformation, we speculate this interval yields a more reliable GPS estimate of 
strain across the Imperial Valley. The GPS rates, as well as those derived 
through conventional geodetic techniques, are listed in Table 5.3. These are 
compared with the predicted relative velocity between the Pacific and North 
American plates (NUVEL-1), and rates derived from VLBI measurements 
between stations along the western coast of California and the stable North 
American continent. The 1988-1989 GPS rate is comparable to the plate 
velocity estimates, while the earlier GPS interval is somewhat higher. This 
suggests that all plate motion is concentrated across the valley, with little or 
no deformation west of the Elsinore fault. Conventional measurements taken 
over the last 50 years indicate significantly smaller rates (3.4-4.3 cm/yr), and 
thus require additional slip on faults not spanned by the networks to satisfy 
the plate velocity. GPS and trilateration (EDM) provide comparable 
accuracies, and both are considerably more precise than triangulation 
(although GPS yields 3-dimensional positions whereas the conventional 
methods do not). Since the EDM observations span a 15 year period 
compared to the 3 year GPS coverage, the trilateration rate should more 
accurately reflect the deformation across the valley, which suggests the GPS 
measurements over-estimate the true displacement. In fact, a 3.4 cm/yr 
deformation rate fits the 1988-1989 observations nearly as well (Figure 5.7). 
Am alternate explanation is accelerated deformation between 1986 and 1989. 
The triangulation data indicate highly time-dependent displacements. 
Between 1941 and 1954 the calculated rate is significantly greater than the 
average between 1941 and 1986, although this is attributed to post-seismic 
effects following the 1940 Imperial Valley earthquake. No increased rate is 
observed following the 1979 earthquake [Savage et ai, 1986J. There is 
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Table 5.3 Displacement Rates 


Method 

Region 

Interval 

Rate 

(cm/yr) 

Reference 

GPS 

Imperial Valley 

1986-1088 

4.8 

This Study 



1088-1980 

5.2 


Triangulation 

Imperial Valley 

1941-1986 

4.3 

Snay and Drew [1988] 



1041-1054 

6.1 




1954-1967 

2.1 




1967-1988 

4.5 


Trilateration 

Imperial Valley 

1972-1987 

3.4 

Prescott tt al. [1987b] 

Plate Model 

Plate boundary 

~ 3 m.y. 

4.7 

DcMcts tt at. (1990] 

VLBI 

Continental 

1970-1987 

4.8-5. 1 

Clark tt al. (1987] 


Kroger tt al. [1987] 
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marginal evidence for a regional strain fluctuation (increase) during 1978 and 
1979 throughout southern California, but the nature of this apparent 
deformation is uncertain [Savage et ai, 1981; Savage et ai, 1986]. Given the 
large uncertainties for the GPS estimates (~ 1 cm), it is not possible with the 

available data to distinguish if there has been increased deformation over the 
last several years. 

Imperial and southern San Andreas fault earthquake potential 

The earthquake recurrence interval along the Imperial fault is estimated 
using the geodetically determined strain rates. The 1940 Imperial Valley 
earthquake ruptured the entire length of the Imperial fault. Approximately 
3.0 and 4.5 m slip (coseismic plus postseismic) are estimated for the northern 
and southern segments of the fault, respectively [Reilinger, 1984]. Surface 
offsets were as great as 6 m south of the border with displacements tapering 
off rapidly to the north [Trifunac and Brune, 1970; Sharp, 1982]. Surface 
rupture was confined to the fault north of the border during the 1979 
earthquake. Geodetic and strong ground motion modeling suggest an average 
slip of about 1 m along the 1979 rupture plane, with patches of higher 

displacement (asperities) [e.g., Hartzell and Heaton, 1983; Archuleta, 1984; 
Reilinger and Larsen, 1986] 

Triangulation measurements along the Imperial fault and north of the 
U.S.-Mexico border indicate concentrated deformation in a narrow 20 km wide 
zone. At an observed strain rate of 4-5 cm/yr and per event ruptures between 
1 and 3 m, a 20-75 year earthquake recurrence interval is calculated for the 
northern Imperial fault. This assumes total strain release during seismic 
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episodes. The estimated recurrence rate is comparable to the 32 year 
earthquake repeat time suggested by Sykes and Nishenko (1984] and the ~ 50 
year interval predicted by Anderson and Bodin (1987). 

Conventional geodetic data north of the Imperial fault indicate 
distributed deformation over a fairly wide region (> 50 km wide) [Snay and 
Drew , 1988]. Presumably the strain that occurs almost exclusively along the 
Imperial fault in the south is transferred by some mechanism to the San 
Andreas, San Jacinto, and Elsinore faults. The slip distribution across each 
fault segment is important from an earthquake hazard standpoint. 

Three Imperial Valley GPS sites have well determined VLBI solutions 
(Very Long Baseline Interferometry) from observations since 1979 ( Clark et ai, 
1987; Sauber et ai, 1989; Ma et al„ 1989). Relative velocities of BLAC-PINY 
and BLAC-MONU for the GPS and VLBI analyses are listed in Table 5.4 and 
illustrated in Figure 5.9. The VLBI velocities vary depending on the 
continental or global nature of each solution and the extent of data 
availability. Only the north-south GPS displacements are considered because 
of the large east-west errors inherent in the 1986 survey. The fault parallel 
velocities are right-lateral assuming relative motions are oriented N40*W. 
The calculated effect from the 1987 Superstition Hills earthquake sequence is 
removed from the GPS data. The VLBI measurements indicate 1.5 to 2.1 
cm/yr fault-parallel (right-lateral) displacement across the San Andreas fault 
(BLAC-PINY) and 3.0 to 3.5 cm/yr across the valley (BLAC-MONU). The 
GPS measurements suggest 1.4 cm/yr displacement across the fault and 3.2 
cm/yr across the valley (Table 5.3 rates are different since they represent 
averages over the entire network). While the BLAC-MONU velocities roughly 
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Table 5.4 Displacement Rates Across San Andreas Fault 


Baseline 

Method 

Interval 

North 

(cm/yr) 

East 

(cm/yr) 

Fault Parallel 
(cm/yr) 

BLAC-PINY 

VLBI 1 

1982-1987 

1.8 

-1.1 

2.1 


VLB I 2 

1979-1988 

1.5 

-1.0 

1.8 


VLBI 3 

1982-1988 

1.2 

-0.9 

1.5 


GPS 

1986-1988 

1.1 


1.4 

BLAC-MONU 

VLBI 1 

1982-1987 

2.3 

-2.7 

3.2 


VLBI 2 

1979-1988 

2.5 

-2.5 

3.5 


VLBI 3 

1982-1988 

2.4 

-1.8 

3.0 


GPS 

1986-1989 

2.5 


3.2 

1 Clark et al. fl987l 


Ma 1988] 
Sauotr [1989] 
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agree with the conventional geodetic measurements (3.4-4.3 cm/yr), the fault- 
crossing displacements are somewhat surprising as they are lower than 
geologic evidence. The long-term geomorphological slip rate along the 
southern San Andreas fault over the last 10,000-30,000 years is estimated 
between 2.3 and 3.5 cm/yr [J Keller et al . , 1982; Weldon and Sick, 1985], with 
2.5 cm/yr a commonly accepted average [e.g., Sieh and Williams , 1990]. The 
geologic slip rate and radiocarbon dating of Holocene offsets along the fault 
suggest a recurrence interval of about 300 years with the last major event in 
1680 [Sieh, 1986]. This logic leads to the conclusion that the potential for a 
major earthquake along the southern San Andreas is high. However, the 
geodetic evidence reported here indicates a smaller strain accumulation rate 
during the last decade, suggesting decreased earthquake potential assuming 
these measurements are indicative of the last few hundred years. This will be 
observed either as a longer recurrence interval or less slip per event. The 
geodetic data are supported by geologic trenching studies which suggest a 
decreasing southern San Andreas slip rate during the past 1000 years Sieh 
[1986]. If this is true, the San Jacinto fault should play a more active role in 
regional tectonics. In fact, the shear strain along the fault determined from 
EDM observations between 1973 and 1984 is nearly the same as that for 
networks which lie on the San Andreas fault Savage et al. [1986]. The two 
fault systems may alternately assume dominant roles in absorbing plate 
motions, as is suggested by variable Quaternary slip rates along the San 
Jacinto fault [SAarp, 1981]. 

Future analyses 

The 1986, 1988, and 1989 Imperial Valley GPS observations are not 
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Figure 5.9: VLBI (solid arrows) and GPS (dashed arrows) velocities at 
stations PINY and MONU relative to station BLAC (Table 5.4). The 
GPS vectors contain large uncertainty in the east-west direction. The 

fault-parallel geodetic velocities across the San Andreas fault are less than 
geologic estimates. 
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sufficiently resolved to accurately map details of the strain distribution in this 
portion of southern California. This is due to poor data quality from the 
1986 measurements, complications due to the 1987 Superstition Hills 
earthquake sequence, the inability to incorporate L2 phase data from several 
sites during 1989, and the short interval spanned by the observations (1 year 
of nonseismic displacement). However, additional GPS data have been 
collected in southern California in March/April 1988 and again in 
February/March 1990 (Figure 5.10). The 1988 measurements were made in 
conjunction with the Riverside County (California) Flood Control District 
and the Riverside County Survey Department. In a network of 62 stations 
spanning an entire cross section of southern California, a maximum of 8 TI- 
4100 GPS receivers were deployed each day over the two week survey. The 
daily observation period was about 4.5 hours. Although the instruments 
collected data every 3 seconds, during the download from receiver to floppy 
disk only every 10th measurement was recorded (30 second epochs). 
Unfortunately, the receivers were not synchronized and the recorded time-tags 
were randomly distributed at 30-second intervals. These data have been 
subsequently processed with the Bernese software. The day to day baseline 
repeatability on the 18 lines with multiple observations is 1.1 cm, suggesting 
that the receiver time-tag ofisets are not a serious source of solution error 
considering the 3 to 5 cm tectonic movements expected in southern California. 

During 1990, a high precision GPS network was established along a 400 km 
segment of the Pacific-North American plate* boundary from the Gulf of 
California in northern Mexico to just south of the junction of the San 
Andreas and San Jacinto faults (~ 34 *N) [R tiling tr tt al., 1990). Twenty- 
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three receiver were reed for upprmcimately two weeks including the TI-4100, 
Trimble 4000SD, and Trimble 4000SDT models. A total of 134 stations were 
occupied during the campaign. Data collection at 103 sites lasted 6 to 7 hours 
each day, while the daily observation interval at 31 stations was 3 to 4 hours 
(half-sessions). The following month (April 1990), 3 additional monuments 
near Upland, California, were established and surveyed by the Riverside 
County Flood Control District in support of university research associated 
with the February 28 (1990) M L 5.5 Upland earthquake. The 1990 
observations are in the process of being analyzed, although obtaining geodetic 
coordinates for the entire survey will be time consuming due to the enormity 
of the data set and the uncertainty in correlating GPS measurements from 

different receivers. 

The 1986-1990 GPS occupation summary for the Imperial Valley, 
Riverside County, and Baja California is listed in Table 5.5. A total of 183 
stations have been occupied at least once (Figure 5.10), and of these 85 have 
been occupied at least twice with most dating back to 1986 or 1988. The 
station coverage does not include kinematic GPS measurements made on 
relatively short transects (few kilometer) crossing the southern San Andreas 
fault [K. Hudnut, personal communication, 1990]. In addition to the dense 
distribution within the Imperial Valley and Baja California, the established 
network north of the Salton Sea provides station coverage extending from the 
California-Arizona border to near the Pacific Ocean. This network will be 
used to constrain the slip distribution along the major fault systems in 
southern California. When combined with geologic data of fault activity, the 
strain estimates will better define the earthquake potential in this region. The 
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Figure 5.10: GPS stations occupied from 1986 to 1990 in Riverside 
County, the Imperial Valley, and northern Baja California. Triangles 
inside circles indicate stations with multiple occupations. The network is 
composed of 183 stations, of which 85 have repeated observations. 



GPS Observations 
1986-1990 
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Table 5.5 GPS Campaign Summary 


Year 

Stations 

Region 

Organization 

1986 

42 

Imperial Valley 

NGS 

1988 

15 

Imperial Valley 

UNAVCO 

1988 

62 

Riverside County 

RCFC /RCSD 

1988 

21 

Imperial Valley 

NGS 

1989 

28 

Imperial Valley 

UNAVCO 

1990 

134 

Imperial Valley /Riverside County 

UNAVCO/RCFC/RCSD/NGS 


NGS - National Geodetic Survey 
UNAVCO - University Navstar Consortium 
RCFC - Riverside County Flood Control 
RCSD - Riverside County Survey District 
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best GPS estimate of secular strain to date is provided by the 11 displacement 
vectors obtain from the 1988 and 1989 surveys. Once the 1990 data are fully 
analyzed, the increased station density (85 stations) and longer measurement 
interval (at least 2 years) will yield an order of magnitude increase in 
deformation extent and resolution. It should be possible to assess the present- 
day strain accumulation rates on the major fault systems in southern 
California to within a few millimeters per year. The network also provides 
good coverage along the southern San Andreas fault, which will be used to 
constrain fault rupture parameters in the likely event of a large earthquake 
within the next several decades [Larsen, 1990]. 

5.6 Conclusions 

GPS measurements from southern California indicate 5.9 ± 1.0 and 5.2 ± 
0.9 cm/yr right-lateral southeast trending displacement across the Imperial 
Valley for the intervals 1986-1988 and 1988-1989, respectively. These rates 
are significantly larger than those obtained from conventional geodetic surveys 
(3. 4-4. 3 cm/yr), suggesting the GPS observations may overestimate the true 
deformation. The earlier measurements contain relatively large errors, and 
are influenced by the 1987 Superstition Hills earthquake sequence. The 1988- 
1989 data are modeled nearly as well by 3.4 cm/yr of valley crossing 
movement. Regardless, a significant secular deformation component is clearly 
observed in the GPS displacements, which is attributed to the relative 
movement between the Pacific and North American plates. There is evidence 
from VLBI and GPS measurements that the strain accumulation rate along 
the southern-most San Andreas fault is smaller than the calculated long-term 
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geoiogic estimate. This indicates a lower earthquake potential for this 
segment of the fault than is presently assumed, and suggests that the San 
Jacinto system plays a more dominant role for relieving the strain 
accumulation in this region. The measurements discussed here are part of a 
larger 183 station GPS network which spans an entire cross section of 
southern California. A total of 134 stations were observed during a recent 
1990 campaign; many of these sites were previously occupied in 1986 and/or 
1988. An order of magnitude increase in resolution and detail regarding the 
strain accumulation rates along the San Andreas, San Jacinto, and Elsinore 
faults is expected once the 1990 GPS data are integrated with the previous 


surveys. 
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GPS Measurements of Deformation Associated with 
the 1087 Superstition Hills Earthquake: 

Evidence for Conjugate Faulting 

Abstract 

Large station displacements observed from Imperial Valley GPS 
campaigns are attributed to the November 24, 1987 Superstition Hills 
earthquake sequence. Thirty sites from a 42 station GPS network established 
in 1986 have been reoccupied during 1988 and/or 1990. Displacements at 
three sites within 3 kilometers of the surface rupture approach 0.5 m. Eight 
additional stations within 20 km of the seismic zone are displaced at least 10 
cm. This is the first occurrence of a large earthquake (Ms 6.6) within a 
preexisting GPS network. Best-fitting uniform slip models of rectangular 
dislocations in an elastic half-space indicate 130 cm right-lateral displacement 
along the northwest trending Superstition Hills fault and 30 cm left-lateral 
offset along the conjugate northeast trending Elmore Ranch fault. The 
geodetic moments are 9.4 x 10 25 dyne-cm and 2.3 x 10 25 dyne-cm for the 
Superstition Hills and Elmore Ranch faults, respectively. Distributed slip 
solutions using Singular Value Decomposition suggest near uniform 
displacement along the Elmore Ranch fault and concentrated slip to the 
northwest and southeast along the Superstition Hills fault. A significant 
component of non-seismic secular displacement is observed across the Imperial 
Valley, which is attributed to interseismic plate-boundary elastic deformation. 
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4.1 Introduction 

The Global Positioning System (GPS) is rapidly becoming one of the most 
important tools to study tectonic deformation. By recording signals from 
earth orbiting satellites it is possible to determine 3-dimensional coordinates 
of geodetic monuments with high accuracy. With repeated observations the 
station displacement or deformation between surveys is measured. GPS can 
be used to monitor secular deformation such as that associated with plate 
motion, or to record rapid strain fluctuations such as those due to seismic and 
volcanic activity. In its final configuration scheduled for the mid 1990’s, 18 
satellites will be deployed in 6 orbital planes (with 3 additional satellites used 
as active spares). When GPS becomes fully operational it will be possible to 
continuously determine 3-dimensional positions anywhere on or near the 
earth. The available satellite constellation existing for the last several years, 
was optimized for North America making geodetic studies in California 
practical. The observation window in which enough satellites have been 
visible to obtain the high accuracies necessary to measure tectonic motion has 
been about 6 to 8 hours each day. 

On November 24, 1987, two large earthquakes separated by 12 hours 
occurred in the northwest portion of the Imperial Valley region of southern 
California. The first event was located on a northeast trending seismic 
lineament and was followed 12 hours later by rupture along the northwest 
trending Superstition hills fault. What makes this earthquake sequence so 
significant from a GPS standpoint, is that it occurred spatially and 
temporally within a preexisting GPS network. This network was established 
in the Imperial Valley in 1986, with partial resurveys in 1988 and 1990. 
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Fifteen stations are located within 20 km of the rupture zone; three stations 

are within 3 km. This is the first occurrence of a large earthquake within a 
preexisting GPS network. 

We compute GPS displace; u the Imperial Valley between 1986 and 
1090. Observed movements of ne*. 0.5 meters are attributed to the 
Superstition Hills earthquake sequence. The earthquake-induced 
displacements are inverted to estimate seismic slip and the corresponding 
geodetic moment along the rupture planes. In addition, there is a large 
component of deformation which can not be explained by the seismic 
disturbance; we assume this to be a manifestation of continuous strain 
accumulation across the Imperial Valley due to the relative motion of the 
Pacific and North American plates. 


4.2 Imperial Valley Seismicity and Tectonics 

The Imperial Valley region of southern California is a complex transition 

tone between crustal spreading in the Gulf of California and right-lateral 

transform motion along the San Andre*, fault (Figure 4.1) | Lomnitz et at, 

1970; Elder, et a A, 1972], The valley is 4-5 million yeare old and has been 

filled by up to 5 km of late Cenozoic sediments (Fats et at., 1982|. The major 

fault systems and structural grain of the valley trend to the northwest, 

roughly parallel to the direction of plate motion. A significant fraction of the 

North American and Pacific relative motion may be accommodated across this 
region. 

The valley is one of the most seismically active regions of California 
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(Figure 4.1) with much of the activity occurring along the Imperial fault and 
within the Brawley Seismic Zone [ Johnson and Hill, 1982). Several large 
earthquakes have occurred in and near the Imperial Valley since 1940 (Figure 
4.1). The Imperial fault ruptured with a M s 7.1 event in 1940 and a M L 6.6 
event in 1979 [U.S. Geol. Surv., 1982). Segments of the San Jacinto fault 
system broke with a Ml 6.2 earthquake in 1954 (Clark segment) and the 1968 
M l 6.5 Borrego Mountain event (Coyote Creek segment). The most GPS 
relevant episode of seismic activity occurred in 1987 along the Superstition 
Hills segment of the San Jacinto fault system, with a M5 6.2 earthquake on a 
northeast trending seismic lineament followed 12 hours later by a M s 6.6 
event on the Superstition Hills fault. 

Conventional geodetic measurements suggest considerable deformation 
across the Imperial Valley. In fact, a significant fraction of the Pacific- North 
American relative plate motion may be accommodated here. New global plate 
model estimates (NUVELrl) [DeMets et al., 1987; DeMets et ai, 1990] predict 
the rate of relative motion between the North American and Pacific plates (at 
Imperial Valley coordinates: 33.0 * N, 115.5 *W) is 4.7 cm/yr oriented 
N39.6 * W. Triangulation measurements spanning this region have been 
modeled as 4.3 cm/yr of plate-boundary deformation averaged between 1941 
and 1986 [Snay and Drew , 1988). Trilateration measurements made by the 
U.S.G.S. between 1972 and 1987 indicate 3.45 cm/yr relative movement 
between stations on opposite sides of the valley [ Prescott et al., 1987b). The 
orientations of the conventional geodetic displacements are approximately 
N40 * W, although to some extent the direction is non-unique and depends 
upon a priori ass um ptions [e.g., Prescott, 1981). In addition, the conventional 
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Figure 4.1: Seismicity and major fault systems of the Imperial Valley. 
The Brawley ' : *mic Zone is the region of anomalously high activity 
between the r Imperial and southern San Andreas faults. The 

valley represents : on zone between crustal spreading in the Gulf 

of California to the : nd right-lateral transform motion along the 

San Andreas system. 
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geodetic measurements indicate that deformation is concentrated in a narrow 
20 km wide zone along the Imperial fault, while to the north it is distributed 
over a region at least 50 km wide. Presumably, deformation is transferred 
from the Imperial fault, which acts as the primary strain release mechanism 
near the border, to distributed shear along the San Andreas, San Jacinto, and 
Elsinore faults. 

4.3 Superstition Hills Earthquake Sequence 

On November 24, 1987 (1:54 GMT), a large M s 6.2 earthquake occurred 
along a northeast trending seismic lineament northeast of the Superstition 
Hills fault (Figure 4.2) [Magistrate et at., 1988]. The focal mechanism and 
aftershock sequence, which extended for 26 km to the northeast and into the 
Brawley Seismic Zone, are consistent with left-lateral strike slip motion on a 
vertical fault. Seven foreshocks were recorded in the 22 minutes prior to the 
main event, including two with M L > 4.0. Surface rupture consisted of a 
complex pattern of left-lateral northeast-trending offsets ranging in length 
from 1.5 to 10 km, and with maximum displacements between 3 and 13 cm 
[Budding and Sharp , 1988; Hudnut et a/., 1989a]. We refer to this northeast 
trending lineament as the Elmore Ranch fault, although more precisely this 
name refers only to the longest of the surface fractures. 

Twelve hours after the Elmore Ranch event (13:15 GMT), a M s 6.6 
earthquake occurred along the northwest trending Superstition Hills fault. 
The epicenter was near the intersection of the Elmore Ranch and Superstition 
Hills faults. Strong ground motion and teleseismic data suggest the rupture 
process for this second event was somewhat complicated, consisting of 
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multiple sub-events [ Bent et al. , 1989; Frankel and Wennerberg , 1989; Hwang 
et al., 1990; Wald tt al., 1990). Surface rupture extended 24 km along the 
previously mapped trace of the fault [Williams and Magistral e, 1989). Up to 
50 cm right-lateral displacement were measured immediately following the 
earthquake. The aftershock pattern was concentrated slightly to the west of 
the fault and did not extend the length of the surface rupture. Magistrals et 
al. [1989] suggested the aftershock sequence was highly correlated with 
basement structure. Both the Superstition Hills and Elmore Ranch events 
triggered sympathetic surface offsets along the Imperial, San Andreas, and 
Coyote Creek faults [McGill et al., 1989; Hudnut and Clark, 1989). 

Significant afterslip was recorded along the Superstition Hills fault 
following the second mainshock. No afterslip was measured along any of the 
surface ruptures associated with the Elmore Ranch event. In fact, all seismic 
activity essentially stopped along this segment after the initiation of the 2nd 
mainshock. 

One of the most interesting aspects of this earthquake sequence is the 
conjugate nature of faulting. That is, two surface ruptures oriented nearly 
perpendicular to each other. As discussed below, this type of fault interaction 
may be typical of Imperial Valley tectonics and may dictate the mode of 
stress/strain transfer from one fault system to another. 

What makes the Superstition Hills earthquake sequence unique from a 
GPS perspective is that it occurred within a preexisting GPS network. 
Stations located near the seismic rupture zone were displaced nearly 0.5 
meters. These movements, as well as smaller displacements observed at 
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Flgure 4.2: Seismicity and surface faulting associated with the November 
24, 1987 Superstition Hills earthquake sequence. A M s 6.2 event occurred 
along a northeast trending structure (referred to here as the Elmore 
Ranch fault) and was followed by 12 hours with a M s 6.6 event along the 
northwest trending Superstition Hills fault. The focal mechanisms, 
aftershock distribution, and surface oflset measurements are consistent 
with left-lateral strike-slip motion along the Elmore Ranch fault and 
right-lateral strike-slip motion along the Superstition Hills fault. A 
significant amount of postseismic slip was observed along the surface 
trace of the Superstition Hills fault, while activity essentially ceased on 
the Elmore Ranch fault after the M s 6.6 event. The shaded strips along 
each fault indicate the geometrical extent of the dislocations used to 
model the geodetic displacements. 






- 174 - 


nearby stations can be inverted to infer properties of the rupture process. 

This is the first time GPS measurements have recorded the deformation from 
a large earthquake. 


4.4 GPS Observations 

After a brief introduction to the Global Positioning System, in this section 
we report how the GPS data were obtained from field campaigns in 1986 , 
1988, and 1990, as well as the processing methods used to obtain station 
coordinates. The GPS displacement vectors between surveys are presented 
and we discuss how measurement errors are handled for this study. More 
complete details about the Global Positioning System, including theoretical 

aspects and processing methods, are found in King et al. [1985], WtlU et al. 
[1987], and Rocken [1988]. 

The signal structure broadcast from each GPS satellite consists of two 
carrier phase signals modulated by a navigational message as well as pseudo 
random codes. The two carrier frequencies, known as the Ll and L2 phases, 
are broadcast at 1575.42 Mhz (Ll) and 1227.60 Mhz (L2). This is equivalent 
to wavelengths of about 19 cm for the Ll and 24 cm for the L2. The 
navigational message contains the satellite coordinates (broadcast ephemeris), 
clock parameters, satellite health, and general system status. The pseudo 
random codes are accurate time marks which allow a GPS receiver to 
determine the transmission time of the signal. When scaled by the speed of 
light, the pseudorange, or the satellitoreceiver distance is computed. If 
measurements to at least 4 GPS satellites are available, and if satellite 
coordinates are known (usually with the broadcast ephemeris), the 3- 
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dimensional receiver position as well as the satellite-receiver time offset can be 
determined. The positioning accuracy with the pseudorange is 1 to 100 m, 
depending on whether the P or C/A code is used, the receiver type, length of 
observation, and the static or kinematic behavior of the instrument. It is the 
pseudorange which will be used for civilian and military navigation. For 
highly accurate geodetic positioning, however, it is necessary to use the carrier 
phase measurements in a post-processing mode. That is, the data collected in 
the field are brought back to the office or laboratory for analysis, usually with 
a fairly robust computer software system. 

GPS Surveys — Data Collection 

The GPS data for this study were collected during 4 Imperial Valley field 
campaigns from 1986 to 1990 (Table 4.1). A total of 46 stations in or near 
the valley have been occupied at least once during this interval (Figure 4.3), 
while 30 sites have been reoccupied since 1986 (Figure 4.3, Table 4.2). TI- 
4100 GPS receivers supporting GESAR software were used during 1986 and 
1988, while Trimble 4000SDT instruments were used during the 1990 survey 
As one of the first commercially available instruments, the TI-4100 is a code 
correlating receiver which allows the simultaneous tracking of up to 4 
satellites at a time. This instrument records the Ll and L2 carrier signals, the 
P and C/A pseudo-codes, and the broadcast message. The Trimble 4000SDT 
is able to record the Ll and L2 phase components for up to 8 satellites at 
once, as well as the C/A code and the broadcast message 

The National Geodetic Survey (NGS) occupied 54 sites in southern 
California in May/June 1986; 42 stations were located in or near the Imperial 
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Figure 4.3a: Imperial Valley and related GPS stations surveyed in 1986, 
1988, and/or 1990. Station names for sites indicated in the shaded box 
are given in Figure 4.3b. 
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Figure 4.3b: Central-southern Imperial Valley and related GPS stations 
surveyed in 1986, 1988, and/or 1990. 
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Table 4.1 Imperial Valley GPS Campaign Summary 


Year 

Month 

Days 

Stations 

Organization 

1986 

Ma> e 

20 

42 

NGS 

1988 

Febr larch 

9 

19 

UNAVCO 

1988 

March /April 

6 

21 

NGS 

1990 

April 

1 

3 

UNAVCO/RCFC 


NGS - National Geodetic Survey 
UNAVCO - University Navstar Consortium 
RCFC - Riverside County Flood Control 
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Table 4.2 Reoccupied GPS Stations 


Station Name 

Abbr. 

1986 

Occupation 
1988 1 1988 2 

Longitude Latitude Elevation 
1990 ' * (m) 

Acute 1934 

ACUT 

• 

• 

• 

■115.6093 33.0300 

-80.91 

A lamrt 

ALAM 

• 

• 

• 

-115.6111 33.1964 

-72.91 

Black Butte NCMN 1982 

BLAC 

• 

• 


• -115.7198 33.6638 

490.01 

Brawley 2 rm 5 

BRAW 

• 


• 

-115.5434 32.9773 

-66.78 

Calexico 1954 

CALE 

• 

• 


-115.5064 32.6645 

- 34.07 

Coach 

COAC 

• 

• 


-115.4070 33.1962 

-8.40 

College 1967 

COLL 

• 

• 


-115.5024 32.8269 

-54.63 

El Centro 2 1959 

elec 

• 


• 

•115.5622 32.7846 

-45.76 

Frink 1934 

FRIN 

• 

• 

• 

-115.6470 33.3603 

-85.39 

GLO Comew 1934 

GLOC 

• 

• 


-115.2465 32.8396 

-15.46 

Hamar 2 1967 

HAMA 

• 

• 

• 

-115.5007 33.0375 

-79.80 

Holt 1924 

HOLT 

• 


• 

-115.3963 32.7814 

-36.79 

Holtville (Alt) 1934 

HLTV 

• 


• 

-115.3821 32.8084 

-38.77 

Imp 1934 

IMPI 

• 


• 

-115.5698 32.8982 

-59.27 

Imperial 1934 

IMPE 

• 


• 

-115.5788 32.8439 

-51.85 

Junction 

JUNC 

• 


• 

-115.0619 32.7092 

8.35 

sj UUv WlUli 

Kane 1939 

KANE 

• 

• 

• 

-115.8237 33.0614 

10.05 

L 589 1967 

L589 

• 

• 

• 

-115.7611 32.9506 

13.62 

Mack 2 1967 bm reset 

MACK 

• 


• 

-115.1441 32.7288 

1.20 

Mello 3 1967 

MELL 

• 


• 

-115.4653 32.7961 

-47.00 

Mound 1934 

MOUN 

• 



• -115.6998 32.9502 

3.40 

Ocotillo NCMN 1982 

OCOT 

• 

• 

• 

-115.7962 32.7901 

-36.43 

Ocotillo 1935 

OCT1 

• 

• 

• 

-116.0017 32.7338 

111.33 

Offset 217 

0217 

• 

• 

• 

-115.3131 32.6782 

-17.16 

Offset 224 

0224 

• 


• 

-115.7055 32.6493 

50.00 

Offset 227 

0227 

• 


• 

-115.8173 32.8414 

70.64 

Orient 1939 

ORE 

• 

• 


-115.4064 32.9168 

-61.36 

Pinyon Flat 

PINY 

• 


• 

-116.4588 33.6093 

1235.88 

T 1226 

T122 

• 


• 

-114.8050 32.8583 

141.19 

Tamarisk 3 1967 

TAMA 

• 

• 


-115.4783 32.8829 

-68.19 


1 UNAVCO 

2 NGS 
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Valley. Most marks were observed at least two days (Table 4.3), although 
redundant baselines were relatively uncommon (i.e., simultaneous occupation 
of the same station-station pair for two or more days). Unfortunately, the 
quality of the observations were very poor. The scheduled 4.5 hour daily 
occupations were somewhat less than the 6-8 hour sessions typical of other 
southern California campaigns. In addition, due to a variety of equipment 
and logistical problems, less than 2 hours of data were collected at several 
sites. In the observation scenario, there was less than a 2 hour period in 
which more than 3 satellites were simultaneously tracked. In fact, there was a 
scheduled 1 hour gap in Satellite 6 during the middle of the measurement 
session. The data were generally noisy and the final double difference phase 
observables may contain uncorrected cycle slips. 

In February/March 1988 university field crews (UNAVCO) observing for 9 
days reoccupied 15 of the Imperial Valley marks, as well as establishing 4 new 
stations along the Salton Sea. These new sites were installed at tide gauge 
monitoring facilities and will be used to constrain vertical GPS baseline 
components and for monitoring tectonic tilting in the Salton Trough. Most 
monuments were occupied for 2 days (Table 4.3) (1 was observed 3 days, 1 for 
4 days, and 2 for 1 day). The scheduled nightly observation scenario lasted 
7.5 hours, with a total of 7 satellites tracked. 

In March/April 1988 the NGS reoccupied 21 of the 1986 stations (7 of 
which were observed by the university crews a month earlier). Most sites were 
occupied only 1 day (Table 4.3). The daily observation period of 6.0 hours 
was slightly shorter than the February /March survey, although 7 satellites 
were tracked each day. 
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Table 4.3 Imperial Valley GPS Occupation History 


Year Day Stations 


1986 141 

COLL 

CALE 

L589 






142 

0217 

0216 

MELL 

E122 

HOLT 

ELEC 

0227 


143 

COLL 

CALE 

TAMA 

HLTV 

MACK 

0224 

IMPE 

0225 

144 

OCTI 

IMPI 

F122 

OCOT 


0224 

0227 

0225 

145 

HAMA 

IMPI 

ALAM 



ACUT 

MOUN 


146 

HAMA 

CALE 

TAMA 

ORIE 


GLOC 

C012 

BRAW 

147 

COAC 

CALI 

ALAM 

ORIE 

FRIN 

BLAC 

KANE 


148 

OCTI 

MONU 

L589 

OCOT 

FRIN 

ACUT 

KANE 

BRAW 

149 

0217 

0216 

J060 

T122 

MACK 

GLOC 

MOUN 

YUMA 

150 

COAC 

MONU 

MELL 

T122 

HOLT 

BLAC 

C012 

YUMA 

151 

G035 

S025 

L589 

D26B 

HOLT 

ACUT 

X35A 

JUNC 

153 

G035 

S025 

J060 

D26B 



X35A 

BRAW 

161 

PINY 


SAND 

OCOT 

PEAR 

BLAC 



162 

PINY 


SAND 

OCOT 

PEAR 

BLAC 


YUMA 

163 

PINY 


SAND 

OCOT 

SANT 

BLAC 


YUMA 

164 

PINY 


JPLP 

OCOT 


BLAC 


YUMA 

167 

PINY 

MONU 

MOJA 

NIGU 

SOLE 



YUMA 

168 

PINY 

MONU 

MOJA 

NIGU 

SOLE 

OTAY 


YUMA 

169 

LAJO 

MONU 

MOJA 

BOUC 


OTAY 


YUMA 

170 

LAJO 

MONU 

MOJA 

BOUC 

CUYA 



YUMA 

1988' 55 

L589 

BLAC 

FTJS 

COAC 

SANl 




56 

ALAM 

GLOC 

OCOT 

SSP1 

SANl 




57 

ALAM 

GLOC 

OCOT 

SSP1 

ORIE 




58 

L589 

GLOC 

BOB1 

COAC 

ORIE 




59 

L589 

BLAC 

FTJS 

KANE 

FRIN 




60 

L589 

CALE 

TAMA 

KANE 

FRIN 




62 

0217 

CALE 

TAMA 






63 

0217 

COLL 

OCTI 






64 

HAMA 

COLL 

OCTI 






1988 2 88 

PINY 

OCOT 

BRAW 

HOLT 

MELL 




89 

PINY 

OCOT 

IMPE 

IMPI 

0217 




90 

PINY 

OCOT 

ACUT 

HAMA 

KANE 




91 

PINY 

OCOT 

ELEC 

HLTV 

OCTI 




92 

PINY 

OCOT 

MACK 

0227 

T122 




93 

PINY 

FRIN 

JUNC 

L589 

0224 




1990 100 

BLAC 

MOUN 

BC01 







Day is Julian day of year 

1 UNAVCO 

2 NGS 
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In addition to the UNA.VCO and NGS campaigns, GPS observations in 
1988 were made at Mojave (California), Westford (Massachusetts), and 
Richmond (Florida). These sites are continuously monitored as part of the 
Cooperative International GPS Network (CIGNET) [CAm, 1988]. Data from 
CIGNET stations are used to improve satellite orbits in the GPS processing. 
Unfortunately, these observations are frequently of poor quality and contain 
abundant cycle slips. During the 1988 campaign, observations were not 
always available or usable at all CIGNET sites on all days (Table 4.4). 

Station MOUN (Mound), which was surveyed in 1986, was destroyed from 
the 1987 earthquake sequence. The site is located less than 1 km from the 
surface rupture of the Superstition Hills fault. Field investigation during 
early 1988 reveiled that the monument and supporting concrete base had. been 
completely uprooted from the ground. Destroyed monuments usually can not 
be tied to previous surveys because of the high accuracy required for crustal 
motion research. If the suspected deformation is significantly large, however, 
some information may be recovered if the monument (or a substitute) is reset 
in approximately the same location. Site inspection at MOUN clearly showed 
where the old monument had been and in early April, 1990 a rebar rod acting 
as a temporary benchmark was set at approximately the same position. We 
estimate the temporary mark was set within 0.15 m of the previous 
monument. Because of its proximity to the 1987 rupture zone, the calculated 
seismic displacement at MOUN is about 0.5 m. Therefore, reoccupation of 

this site should retrieve a tectonic component larger than the expected 
uncertainty. 

The 1990 survey was conducted to establish the displacement of MOUN 



- 186 - 


relative to its 1986 position. This "mini-campaign" (Table 4.3), which 
included stationary GPS receivers at- only three sites (MOUN, BLAG, and 
BC01) was performed interactively with kinematic GPS transects along the 
southern San Andreas fault (IC. Hudnut, personal communication, 1990). 
Data were collected for only one night; a total of 9 satellites were visible 
during the scheduled 7.0 hour experiment. 

GPS Processing 

The 1986 GPS observations were processed with the GPS22 software 
developed at the National Geodetic Survey. Satellite orbital information was 
provided by the NSWC (Naval Surface Weapons Center). A tropospheric 
delay parameter was solved at each station, constrained by surface 
meteorological measurements. Ambiguities were fixed to the nearest integer. 
Each of the 20 days of observation was processed separately, and the daily 
solutions were combined to form one set of station coordinates with the 
geodetic adjustment program DYNAP (DYNamic Adjustment Program) 
[Drew and Snay, 1989). All coordinates were computed in the WGS-84 
reference frame [Defense Mapping Agency, 1987). Due to poor data quality 
and observation constraints, the accuracy of these measurements is believed to 
be on the order of 1 ppm (parts per million) [ Neugebauer , 1988). 

Data from both 1988 campaigns were processed with the Bernese GPS 
analysis software (version 3.0), from the University of Bern in Switzerland. 
(Preliminary processing such as data translation and cycle slip fixing were 
performed with an earlier version of the software.) The Bernese GPS analysis 
package allows 3-dimensional station coordinates to be determined from the 
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integration of carrier phase, pseudo-range, and orbital data [e.g., Gartner et 
al., 1985; Bevtlcr et al ., 1985; Rocken, 1988]. In addition to measuring station 
coordinates, the double difference estimation algorithm can solve for 
adjustments to six Keplerian elements and two radiation pressure coefficients 
for each satellite, an atmospheric zenith delay parameter at each station, a 
clock error term at each site, ionospheric model coefficients, and cycle 
ambiguity terms [e.g., Rocken , 1988]. 

For each of the 1988 surveys, all data were combined into a simultaneous 
multi-day solution. Surface meteorologic data (temperature, pressure, relative 
humidity) were used to constrain a Saasmotonian atmospheric model, and 
independent tropospheric zenith delay parameters were estimated at each 
station. We experimented with fixing ambiguities but found mixed results, 
therefore, ambiguities were left unresolved in the final solutions. 

In addition to station coordinates, satellite orbital parameters were 
estimated for both solutions. GPS observations from the CIGNET tracking 
sites were used to constrain the orbits. The Bernese software is able to 
combine ephemeris information from several days into a single multi-day arc. 
It has been suggested that multi-day satellite arcs sig nifican tly improve GPS 
precision [e.g., Ltchten, 1987]. Typically, 6 Keplerian elements and 2 
radiation pressure coefficients are estimated for each satellite. Davis et al. 
[1989] utilized 4-day arcs (72 hours plus the length of the daily observation 
session) to analyze GPS data from North America. A 5-day satellite arc was 
used to process GPS data off the southern California coast near Santa 
Barbara [ Larsen et al ., 1990], 
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There comes a point, however, when the force model describing the 
satellite orbits begins to break down. It is unlikely that the 10-day interval 
of the February/March 1988 survey (9 observation days plus one day of no 
measurements) could be modeled from a 10-day satellite arc. For purposes of 
orbit determination, therefore, this campaign is divided into multiple 3-day 
arc segments (days 55-57, 58-60, and 62-64), each defined by its own initial 
conditions. The data are still processed simultaneously, except that 
independent satellite parameters are estimated for each 3-day arc. For the 
entire campaign, 18 Keplerian orbital elements and six radiation pressure 
terms are estimated for each satellite. A similar technique is used for the 
March/ April survey. The satellite orbits for this campaign are defined by two 
3-day arcs (days 88-90 and 91-93). 

Orbits are improved by holding the coordinates of at least 3 CIGNET 
stations fixed, to values well determined from VLBI and satellite laser ranging 
[Murray and King, MIT interoffice memorandum, 1988]. The GPS phase 
observables from these fiducial sites constrain the orbits, which in turn 
improve the solution accuracy of the unfixed stations estimated as free 
parameters. Notice in Table 4.4 that there are several gaps in the fiducial 
coverage. Although it would be better to have complete data, in general, this 
is not a problem since there is always at least one day in each 3-day arc 
interval where observations from 3 fiducial sites are available. 

After each 1988 campaign was processed separately, the Cartesian 
coordinate differences from the two solutions were adjusted by least squares to 
obtain one set of coordinates for that year. 
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Table 4.4 CIGNET Data Availability 1988 
Day MOJA WEST RICH 



Day is Julian day of year 
• Good quality data 
□ Poor quality data 
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The 1990 "mini- campaign" was processed with the Trimvec software, 
made available from the receiver manufacturer (Trimble Navigation, 
Sunnyvale, California). Recall that the purpose of this survey was to 
establish the displacement at station MOUN, presumed large because of its 
proximity to the 1987 seismic rupture. Since this mark had been reset 
between the 1986 and 1990 occupations, the error in the displacement 
estimate will be fairly large (~ 15 cm). Therefore, high accuracy from a GPS 
perspective is not required. Only stations BLAC and MOUN were included in 
the processing since BC01 was not occupied in 1986. The orbits were given 
by the broadcast ephemeris, and surface meteorological data were used to 
constrain a tropospheric delay model. (This line was processed by K. 
Hudnut.) 


Station Displacements 1080-1088 

GPS station displacements for the interval 1986 to 1988 are shown in 
Figure 4.4. All movements are made relative to station OCTI. Only 
horizontal components are shown. The method in which errors are 
formulated and utilized is disussed below. Generally, the observed 
displacements can be decomposed into 3 components: l) seismic deformation 
due to the Superstition Hills earthquake sequence, 2) secular deformation due 
to the Pacific-North American relative plate motion, and 3) measurement 
error attributed to poor data quality from the 1986 survey, most notable in 
the east-west direction. 

The GPS displacement vectors suggest considerable deformation between 
the 1986 and 1988 campaigns, of which a significant fraction is attributed to 
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Figure 4.4: Imperial Valley GPS station displacements between 1986 and 
1988. All movements are relative to station OCTI. The observed 
displacements are attributed to the 1987 earthquake sequence, secular 
pi ate- bound ary deformation across the Imperial Valley, and measurement 
error. Movements near the 1987 rupture zone approach 0.5 m. Error 
ellipses are determined by multiplying the formal errors by a variance 
factor, determined so the average error scales as 1 ppm (parts per 
million). The uncertainty in the east-west direction is about 4 time larger 
than that in the north-south direction. The anomalous southwest 
trending apparent movements for those stations to the southeast are 
attributed to measurement error. 
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the 1987 Superstition Hills earthquake sequence. Stations nearest the seismic 
rupture zone (KANE and L589) show movements on the order of 40 cm. In 
fact, the 13 km KANE-L589 baseline was shortened by 70 cm. The 
orientations of the displacements are consistent with the conjugate fault 
pattern indicated by the mapped surface offsets (i.e., right-lateral rupture on 
the Superstition Hills fault and left-lateral rupture along the Elmore Ranch 
fault). Other stations near the active fault system appear to have been 
affected by the 1987 event as well. 

There is an additional component of displacement not readily explained 
by the seismic deformation. Stations east of the Imperial fault tend to be 
moving south or southeast relative to sites on the other side of the valley. 
This secular displacement is attributed to the relative motion between the 
North American and Pacific plates. Unfortunately, it is difficult to ascertain 
the magnitude of this deformation since station coverage west of the Imperial 
fault is somewhat lacking and many of these sites have rather large 
seismically induced displacements (recall that all GPS vectors are relative to 
OCTI). However, the pi ate- boundary deformation does appear to be 
considerable, which is not too surprising considering conventional geodesy 
indicates a significant fraction of relative plate motion is occurring across the 
Imperial Valley. 

Also evident in Figure 4.4 are unusual movements which do not appear to 
be tectonically related. Most notable are the southwest trending vectors (as 
opposed to southeast) for those sites near the border east of the Imperial 
fault. It appears as if the entire network has undergone a systematic- 
clockwise rotation. We have investigated this possibility by assuming the 
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network could be rotated (and translated) in terms of an outer coordinate 
solution by minimizing the displacement component perpendicular to the 
structural axis of the valley (N40*W) [e.g., Prescott, 1981]. Stations KANE 
and L589 were not included in the solution. However, the applied adjustment 
did not correct for the anomalous displacements, and in fact, made the 
apparent deformation less uniform. Although these unusual movements can 
not be attributed to a simple coordinate rotation, they can be explained by 
large east-west trending systematic measurement errors in the 1986 survey. 
This is consistent with the longitudinal orientation of the computed error 
ellipses (see below) and suggests the north-south displacement components 
may be a more reliable indicator of tectonic deformation. 

Station Displacements 1980-1090 

The 1986-1990 displacement of MOUN relative to BLAC is shown in 
Figure 4.5 (dashed arrow). Although the errors are large because MOUN was 
reset between surveys, the GPS data indicate significant movement attributed 
to the Superstition Hills events (recall MOUN is less than 1 km from the 1987 
surface rupture). However, based on conventional geodesy as well as the 
1986-1988 GPS movements, a non-seismic plate-boundary displacement 
component is suspected in the measurements. We attempt to remove this 
component by estimating the MOUN-BLAC secular displacement based on 
EDM observations. This is discussed in more detail below. We remove 2 
years of the accumulated MOUN-BLAC secular movement, and then compute 
the displacement relative to station OCTI. Although the measurements still 
contain 2 years of non-seismic deformation, it places the displacements into 
the same reference frame as the 1986-1988 movements. The adjusted MOUN 
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Figure 4.5: Imperial Valley GPS station displacements between 1986 and 
1990. The displacement at MOUN relative to BLAC is shown by the 
dashed arrow. For consistency with the 1986-1988 observations the 
estimated displacement at MOUN relative to OCTI is calculated by 
subtracting the estimated secular velocity of BLAC relative to OCTI 
obtained from conventional geodetic measurements (see Figure 4.9). 
Station MOUN was reset between surveys and site inspection during 1990 
revealed a positioning error of about 15 cm. 
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displacement is shown in Figure 4.5 (solid arrow). 

Station Displacements 1088-1088 (February- April) 

Seven GPS sites were occupied during both 1088 campaigns (Table 4.2). 
Calculated site displacements for this 1 month interval are shown in Figure 
4.6. An adjustment (simple translation) was applied to all movements so that 
the sum of the vector displacements is zero. The magnitude of the apparent 
movements range from 0.9 to 2.9 cm, averaging 1.6 cm. It is interesting that 
the 2.9 cm displacement at KANE is to the west-southwest, more or less 
expected if left-lateral afterslip occurred along the Elmore Ranch fault. 
However, while postseismic offsets were significant along the Superstition Hills 
fault, almost all Elmore Ranch activity ceased after the initiation of the 2nd 
main event [Williams and Magistrals, 1989; Magistrals, 1989; Hudnut st al. , 
1989a]. The observed displacements probably represent measurement error as 
opposed to real deformation. Note that the largest vector components are 
oriented in the east-west direction. Station 0217 also exhibits a fairly large 
apparent east-west directed movement, although this is almost certainly not 
tectonically related. 

Even if the displacements shown in Figure 4.6 are entirely measurement 
error, they do illustrate two points: 1) the accuracy with GPS is easily 
sufficient to monitor tectonic motions, and 2) the anomalous displacements for 
the 1986-1988 interval (Figure 4.4) are probably due to the poor quality of the 
1986 data. Typical crustal deformation rates across major tectonic structures 
in southern California are between 1 and 5 cm/yr. At 1 to 2 cm GPS 
accuracy, such deformation would be resolvable in time scales as short as 1 
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year. Although the short baseline (~ 50 km) horizontal precisions computed 
from GPS repeatability studies are generally at the sub-centimeter level [e.g., 
Dong and Bock, 1989], these tests usually involve multiple occupations of the 
same network over a consecutive 4 to 5 day interval. Because the 7 stations 
shown in Figure 4.6 were not all observed simultaneously (Table 4.3), the 
repeatability is somewhat degraded since errors will tend to propagate 
through the solutions. However, the relatively good consistency suggested 
with the 1988 results suggest that some of the unusual movements observed 
with the 1986-1988 displacement vectors are probably due to poor data 
quality from the 1986 campaign. 


GPS Errors 

Formal estimates of GPS uncertainty almost always underestimate 
variances derived from repeatability studies. We attempt to determine more 
realistic and illustrative errors by multiplying the structure of the formal 
covariance matrix calculated with the GPS solution by an estimated variance 
factor, which scales as average baseline length. For the 1986-1988 
displacements, the primary error source is presumably due to the 1986 survey, 
where the estimated accuracies are on the order of 1 ppm. We have chosen a 
variance factor so that the average baseline error scales as 1 ppm. Since the 
1988 data were processed robustly utilizing orbit improvement techniques (not 
available with the 1986 data), we assume these errors are negligible. 
Although this method is somewhat ad hoc, it does illustrate an important 
fundamental. Notably, the error in the east-west component is about 4 times 
larger than the north-south uncertainty. This is attributed to the north- 
south ground-track of the satellite orbits. 
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Figure 4 . 6 : Imperial Valley GPS station displacements between 
February /March 1988 and March/April 1988. The observed movements 
indicate the magnitude of errors due to the 1988 survey. The vector scale 
is twice that of Figures 4.4 and 4.5. The displacement at KANE could 
represent postseismic deformation from the 1987 earthquake sequence. 
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Be cause station MOUN was reset between the 1086 and 1990 surveys, the 
largest error source is due to the difficulty in establishing the new mark at the 
previous location. From field inspection, this uncertainty is estimated at 15 
cm (in all directions) and is additive to the GPS measurement error. 

4.5 Modeling 
Theory 

Simple dislocation theory is often used to model seismically induced 
geodetic deformation. The earth is considered a homogeneous isotropic elastic 
half-space with no stress applied to the free surface. The displacement field 
u k for a dislocation £ in the medium is given by 

uk- i^ vu ‘ w fr dE (4.i) 

where Vu, is the discontinuity, w* are the displacement Green’s functions due 
to a set of strain nuclei, and Vj are the direction cosines of the normal to the 
surface element d£ [ Steketee , 1958; Chinnery, 1961]. Analytical solutions to 
this integral are rather complex, but have been simplified for special cases of 
dislocation or fault geometry (e.g., Chinnery , 1961; Savage and Hastie, 1966; 
Mansmha and Smyhe, 1967). General expressions of the displacement field for 
rectangular strike and dip-slip faults of arbitrary inclination have been 
computed by Maneinha and Smylie [1971] and Okada [1985]. Arbitrary slip 
directions can be designed by the superposition of strike and dip-slip 
dislocations. 

The strain/stress within a medium is computed by differentiating the 
displacement field. For the displacement u k where u is a function of the 
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geometrical coordinates x„ the components of the strain tensor E are given by 

+ (4.2) 

11 2 [dxj d*i j 

In an isotropic medium the stress tensor <r is given by 

2^j ( 4 * 3 ) 

where 6 is the dilatation (*-£«„). We co“s ider the medium 83 a Poisson 


solid with \-m-2.8x 10 11 dyne-cm 2 . 


A rectangular dislocation within an elastic half-space will create a 
spatially dependent stress tensor a throughout the volume. The force acting 
at a point along an arbitrarily oriented plane in the medium is computed by 
multiplying cr by the outward normal vector to the plane (N n ). That is, the 
traction vector T is given by 

T-ofl n (4-4) 

If we assume the plane is coincident with a fault, then the forces generated on 

this secondary structure due to the initial dislocation are determined by 
calculating the traction vectors at selected points along the fault. The normal 
(<r n ), strike-slip (<7 S ), and dip-slip (tr d ) stresses on the fault plane are computed 

by 


<r n -T-fl n 

(4.5a) 

ct s -T-K 8 

(4.5b) 

<r d -T*fl d 

(4.5c) 


where fl B , H s , and K d are the normalized vectors perpendicular to, along 
strike, and along dip to the fault plane. Analytic solutions for the dislocation 
generated stress and strain fields within a medium are given by Iwasaki and 
Sato (1979] and Altwinc [1974]. 
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Inverse Methods 

Inversion of seismically generated geodetic displacements can yield fault- 
rupture parameters, such as the slip distribution along a fault plane [e.g., 
Ward and Bamentos, 1386; Harris and Segall, 1987; Segall and Harris , 1987; 
Snay, 1989). We use a nr od similar to that outlined in Segall and Harris 
[1987]. Singular Value De osition (SVD) [e.g., Lanczos , 1961; Jackson, 
1972; Menke , 1984) and elastic dislocation theory are used to invert the 
Imperial Valley GPS measurements for seismic slip along the Superstition 
Hills and Elmore Ranch faults. 

The relationship between surface deformation and slip along a rectangular 
dislocation is defined by Equation (4.1). The rupture plane is modeled as a 
set of non-overiapping rectangular dislocations. That is, the fault plane is 
partitioned into multiple sub-eiements or sub-faults. The slip distribution 
along the seismically active fault is given as the discreate approximation of 
slip along each sub-element. The normal equations which govern surface 
displacement resulting from such slip is given by 

A«m f = d* (4. 6 ) 

where the superscripts g and / refer to geodetic observation and fault slip, 

respectively. Each row of A 8 is determined from (4.1), and is a function of 

sub-fault geometry and geodetic position. The slip distribution m f is defined 

by m f *»[mj,m 2 , ■ • • ,m,J T , where m 4 is the slip along the tth sub-fault. The 

data vector d 8 contains the geodetic observables. Each GPS station 

displacement will add 3 rows to A 8 and 3 elements to d 8 , corresponding to the 

vertical and two horizontal components. In practice, we may choose to ignore 

the less accurate vertical observation. This is especially true in strike-slip 



- 203 - 


environments where the predominant displacement direction is horizontal. 


Surface rupture is easily included into (4.6) by considering measurements 
of surface displacement as geodetic observation. The offsets are modeled as a 
priori slip information on the surface intersecting sub-faults. Equation (4.6) 
becomes 



where d, 3 are the discreate approximations of surface slip along the fault trace 
and A, 3 = 1 if sub-fault element j corresponds to surface slip offset t; otherwise 
A*=0. 


The GPS displacements shown in Figures 4.4 and 4.5 are not connected to 
an external reference but are defined relative to station OCTI. (Some geodetic 
measurements such as line-length change recorded by EDM observations are 
independent of an absolute reference.) The displacement at OCTI is assumed 
to be 0, which in fact may be true from a seismic standpoint since this site is 
far from the Superstition Hills rupture zone. However, any attempt to use 
the GPS displacements as a criteria for evaluating the effect of the earthquake 
sequence will be distorted by measurement error at OCTI and/or non-seismic 
deformation between this site and the other stations. This ambiguity is 
largely circumvented if displacement-offset terms are estimated in addition to 
the fault slip parameters. Equation (4.7) is then rewritten 


Am “A 


m* 

m° 


(4.8) 


where m,° is the ith non-seismic component (i.e., north-south, east-west, 
vertical) uniformly added to all station displacements. 


The Singular Value Decomposition of A is given by 

A-UXV t ( 4 . 9 ) 

where U is a matrix of eigenvectors spanning the data space, V is a matrix of 

eigenvectors spanning the parameter space, and X is a diagonal matrix of 

singular values. Without loss of generality this is written 

■^.“tJpXpVjT (4.10) 

where p refers to the non-zero singular values. If the normal equations of 

(4.8) are normalized to have unit variance [e.g., StgaU and Harris, 1987), the 

generalized inverse of (4.8) and (4.10) is given by 

A " 1==V pV 1U p T (4.11) 

[Lanczos, 1961; Menkc, 1984). In practice it is often necessary to restrict the 

volume of the parameter space by considering only the k largest singular 

values, setting all others to 0. 

The generalized solution to (4.8) for the k largest singular values is given 
by 

m-Afc-M+VoOb (4.1,2) 

where V 0 are eigenvectors spanning the null space of the model and Qq is a 

vector of arbitrary coefficients. The volume of the model space not 

constrained by observation is defined by V 0 o^. This term is not influenced by 

the geodetic data and is thus arbitrary. Often it is the minimum-length 

solution m — Af *d which is of interest (the coefficients of Oq are 0). However, 

some other solution criteria can be satisfied by carefully designating the 

coefficients of Qq. 

For high resolution fault models where the rupture plane is partitioned 
into numerous sub-faults, it is necessary to apply some type of smoothing 
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constraint over the dislocation surface to prevent the slip distribution from 
taking on an oscillatory pattern. Segall and Harris [1987] showed that the 
"roughness" of fault slip could be minimized by considering smoothness as an 
a priori constraint utilized from the model null space through the coefficients 
of Oq (Equation 4.12). They considered a smoothing matrix T, with 
coefficients determined from the discreate approximation of the Laplacian 
operator V^c^/dx 2 +<9 2 /d y 2 , where x and y are the fault distances along 
strike and dip, respectively. The boundary conditions around the lower and 
lateral edges of the dislocation are assumed to be null slip, so that the applied 
smoothing operator causes the calculated fault offset to tend toward zero 
along these boundarys. Because the Superstition Hills and Elmore Ranch 
faults ruptured the surface, the upper boundary is considered an 
unconstrained dislocation. The estimated fault slip is then given by 

m - [I - V 0 (V 0 T T T TV 0 )- 1 V 0 T T T T]A k - 1 d (4.13) 

(Equation 13, [ SegaU and Harris , 1987]). A similar formulation for utilizing 

fault smoothness over the model null space is given by Harris and Segall 

[1987]; an alternate method considering fault smoothness as quasi-data is 

provided by Snay [1989]. 

For over-constrained solutions, where there are more independent data 
than parameters, if k = p then SVD is equivalent to simple least-squares. 
This is advantageous since the solution provided by (4.12) can be utilized for 
either uniform dislocations or for detailed parameterizations where the fault 
plane is partitioned into multiple sub-elements. 

Simple dislocation theory has the advantage that the displacement and 
stress/strain fields for simple fault ruptures can be computed almost 
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instantaneously. The inverse problem of using geodetic data to calculate the 
slip distribution along a fault plane is also straightforward. However, we 
have assumed the earth can be modeled as a homogeneous half-space. Crustal 
layering or inhomogeneities in the earth can introduce non-existent structure 
into half-space models [Savage, 1987). While low- resolution schemes such as 
the average slip over the fault plane will not be seriously affected, attempts to 
resolve detailed properties may be badly contaminated by artifacts of earth 
structure. 


Seismic Displacement 

We model rupture along the Superstition Hills and Elmore Ranch faults 
as strike-slip dislocations along vertical planes extending from the surface to 
10 km depth. Each dislocation approximately coincides with the mapped 
surface rupture and/or aftershock distribution. The geometrical parameters 
for the modeled faults are listed in Table 4.5. Initially, the Superstition Hills 
and Elmore Ranch faults are considered uniform dislocations, with no slip 
variation allowed on the rupture planes. The dislocations are then 
partitioned into multiple sub-elements and the slip distribution along the two 
faults is calculated from the discrete offset for each sub-fault. 

For an initial estimate of slip along the two faults, we consider only the 
displacements at KANE, MOUN, and L589; the three GPS sites nearest the 
seismic rupture zone (Figure 4.7). Because the observed movements at these 
stations are relative to OCTI, we also solve for a uniform north-south and 
east-west offset in the displacements. This will remove any systematic 
distortion due to measurement error at OCTI. Both faults are regarded as 
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Table 4.5 Parameters for Modeled Dislocations 



Superstition Hills 

Elmore Ranch 

Length (km) 

25 

25 

Width (km) 

10 

10 

Strike 

N50* W 

N40 * E 

Dip 

90 

90 

Depth (km) 

0 

0 

Latitude ( * N) 

32.9569 

33.1078 

Longitude ( * E) 

-115.7431 

-115.7505 


Latitude and Longitude are coordinates 
at top center of dislocation 
Depth is depth to top of fault 
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Figure 4.7: The best-fit solution to Model 1 (3-station inversion). The 
solid arrows indicate the observed displacements, while the dashed arrows 
represent the computed displacement based on Model 1. The shaded 
region indicates where horizontal displacements are greater than 4 cm. 
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simple dislocations not segmented into sub-regions. The solution to this 
simple model is given in Table 4.6 (Model 1) and the observed vs. calculated 
displacements are shown in Figure 4.7. Generally, there is good agreement 
between model and observation. This is not altogether unexpected since we 
are solving for 4 parameters with 6 data. The large discrepancy at MOUN is 
presumably due to the uncertainty in relocating the 1990 reset monument at 
the 1986 position (recall the estimated error is 15 cm in each direction). The 
surface deformation pattern computed from Model I is fairly extensive; the 4 
cm horizontal displacement contour is shown in Figure 4.7. 

The north-south and east-west displacement components at each Imperial 
Valley site, plotted as a function of distance from a N40*W trending line 
through OCTI are shown in Figure 4.8. Shown are the non-seismic 
movements; that is, the seismic component computed from Model 1 is 
subtracted from the observed displacements (Figure 4.4). Stations where the 
seismic correction is greater than 4 cm (in the component plotted) are shown 
as open circles, other sites as filled circles. The displacements represent a 
cross-section of non-seismic deformation perpendicular to the plate motion 
direction. 

A fairly consistent pattern is observed in the north-south components 
(Figure 4.8). Stations to the northeast have moved south (or southeast) 
about 8.1 cm relative to sites on the other side of the valley. Stations which 
display the largest deviation are for the most part those sites where the 
applied seismic correction is greater than 4 cm (open circles). This may 
indicate additional fault complexity not accounted for by the simple 
dislocation model used to remove the effects of the 1987 earthquake. The 
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Table 4.6 Inverse Models of Seismic Slip 


Model 

Fault 

Sub-faults 
Strike Dip 

Slip 

(cm) 

Moment 
(xlO 25 dyne-cm) 

Model 1 

SH 

1 

1 

109. ± 13. 

7.8 

Model 1 

ER 

1 

1 

-45. ± 19. 

3.4 

Model 2 

SH 

1 

1 

130. ± 8. 

9.4 

Model 2 

ER 

1 

1 

-30. ± 10. 

2.3 

Model 3 a 

SH 

10 

5 


9.9 

Model 3a 

ER 

10 

5 


5.9 

Model 3b 

SH 

10 

5 


8.4 

Model 3b 

ER 

10 

5 


7.0 

Model 3c 

SH 

10 

5 


6.2 

Model 3c 

ER 

10 

5 


3.9 

Model 3d 

SH 

10 

5 


9.2 

Model 3d 

ER 

10 

5 


4.9 


SH - Superstition Hills fault 
ER - Elmore Ranch fault 
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Figure 4.8: The north-south and east-west GPS displacement 

components for the 1986-1988 interval. All distances are computed 
relative to OCTI on a cross section trending N50“E, perpendicular to the 
plate motion direction. The effects of the 1987 Superstition Hills 
earthquake sequence based on Model 1 are removed. Open circles indicate 
stations where the seismic correction is greater than 4 cm (for that 
component). The 8.1 cm north-south offset between stations on opposite 
sides of the valley is equivalent to 5.9 cm/yr displacement oriented 
N40 * W. The large scatter for the east-west component is presumably 
due to large measurement errors in the 1986 survey. 
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east-west oriented displacements, however, show large data scatter; no 
distinguishable pattern is readily visible across the valley. The scatter is not 
a function of the magnitude of the applied seismic correction, so presumably 
it represents the large measurement errors inherent in the east-west direction. 

We assume the 8.1 cm north-south offset (Figure 4.8) is attributed to 
plate-boundary deformation tween the North American and Pacific plates. 
Taking into account the c. ^ * of the suspected deformation (N40’W), 
as well as the time interval beiw^u „ae 1986 and 1988 survey* 5 years), 
the north-south movements are consistent with 5.9 cm/yr displacement across 
the valley. This is significantly larger than the 3.45-4.3 cm/yr rates obtained 
from conventional surveys. Although accelerated deformation between the 
GPS campaigns can not be ruled out, there is relatively poor station coverage 
in the southwest portion of the valley so it is difficult to estimate the valley 
crossing displacement precisely. TL'.s is even more true considering most of 
the southwestern sites suffered large seismic displacements during 1987. If 

Model 1 (Table 4.6) does v.o* ^ V reflect the rupture process during the 

Superstition Hills earthq the unmodeled effects will propagate 

into the non-seismic estimate. 

The fault rupture calculated from Model 1 depends heavily on shallow slip 
since the three stations used are all within close proximity to the dislocation 
planes. While the observed surface offsets indicate rupture extends to the 
surface, greater slip at depth may go undetected. Therefore, it is necessary to 
examine the displacement at stations away from the fault to ascertain the 
depth extent of faulting. 
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To incorporate more GPS data into the fault-plane inversion, it is 
necessary to remove the non-seismic deformation from the displacement field. 
Perhaps the best example of secular deformation across the Imperial Valley is 
provided by U.S.G.S. trilateration measurements between 1972 and 1987 [e.g., 
Prescott ct al. , 1987a; Prescott et ai, 1987b]. Computed station velocities for 
the Sal ton Trough EDM network (Figure 4.9) are roughly parallel to the 
direction of plate motion (N40'W), although to some degree the geodetic 
orientation is dictated by the outer coordinate solution imposed to transform 
EDM line-length changes into station displacements [ Prescott , 1981], The 
total differential velocity across the network is 3.45 cm/yr and is 
accommodated in a 50 km wide zone [ Prescott et al., 1987b]. Triangulation 
measurements suggest a larger rate between 1941-1986 (4.3 cm/yr) but these 
observations are not as accurate as the EDM measurements [Snay and Drew, 
1988]. However, along the Imperial fault where EDM sites are sparse, the 
triangulation data suggest concentrated deformation in a narrow 20 km wide 
zone. The conventional geodetic data are modeled using the following 
empirical approach. The differential velocity across the valley is taken to be 
3.45 cm/yr. Running along the axis of the valley is a transition zone(s) 
(Figure 4.9), where the strain gradient is defined by simple shear with the 
displacements oriented N40*W. The transition zone north of the Imperial 
fault is 50 km wide; to the south it is 20 km wide. The modeled station 
movements shown by the dashed arrows in Figure 4.9 fit the observed EDM 
displacements extremely well. 

The secular deformation derived from the conventional measurements is 
removed from the observed GPS displacements (Figures 4.4 and 4.5) leaving 


- 210 - 


Figure 4.0: Imperial Valley EDM station velocities computed between 
1972 and 1987 (solid arrows). The movements are largely attributed to 
secular deformation due to the relative motion between the North 
American and Pacific plates. The displacements are modeled (dashed 
arrows) by 3.45 cm/yr displacement across the valley, with right-lateral 
simple shear oriented N40 * W occurring in a transition zones 50 km wide 
north of the Imperial Fault and 20 km wide to the south (shaded). The 
secular deformation is subtracted from the GPS displacements shown in 
Figures 4.4 and 4.5. 
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the seismic component (and measurement error). For Model 2 (Table 4.6) the 
uniform slip along the Superstition Mils and Elmore Ranch faults is 
recomputed using all GPS data (without the secular deformation). The 
residuals (observed minus calculated) for this model are shown in Figure 4.10. 
The largest station discrepancies between model and observation trend in the 
longitudinal direction and are especially noticeable for those sites in the 
southeast. This simply reconfirms our speculation for large east-west trending 
errors. However, the residuals at the three stations nearest the seismogenic 
tone are unusually large. The large vector at MOUN is easily explained since 
this station was reset between surveys. The residuals at LS89 and KANE, 
however, are significantly larger than the average discrepancy computed for 
the other stations. Because both sites are located in close proximity to the 
earthquake rupture tone, this suggests additional seismic slip not accounted 
for by the simple dislocation parameters used for Model 2. 

Although the inversion results (Table 4.6) between Model 1 and 2 are 
marginally different (less slip on the Elmore Ranch fault; more on the 
Superstition Mils fault), this is not significant considering the estimated 
uncertainties. In fact, because the near-field (Model 1) and far-field (Model 2) 
solutions are similar, this suggests that to first order there is not significant 
Slip dependence with depth (within a factor of 2 or 3). It is also noteworthy 
that the uncertainties improve by only ~ 50 % with the additional data 
supplied with Model 2. This illustrates the uec«sity of measurement, near 
the seismic rupture. For simple fault models, where uniform slip is 
constrained to a single dislocation plane (or two plane, in the case of the 
Superstition Mils sequence), it is as important to have at least m in im .i 
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station coverage within a few kilometers of the seismogenic zone as it is to 
have many sites located away from the fault(s). 

Seismic Slip Distribution 

Uniform dislocations along the Superstition Hills and Elmore Ranch faults 
were utilized for Models 1 and 2. To estimate the seismic slip distribution, it 
is necessary to partition the rupture planes into multiple regions or sub-faults 
(Model 3). The divisions must be sufficiently dense as to provide reasonable 
slip resolution. We choose 10 sub-fault elements in the horizontal and 5 in 
the vertical, so that each fault is partitioned into 50 sub-regions. The 
dimensions of each dislocation element is 2.0 km in width (vertical) and 2.4- 
2.6 km in length (2.4 km for Superstition Hills and 2.6 km for the Elmore 
Ranch fault). The slip distribution along the two rupture planes is 
constrained to be sufficiently smooth (Equation 4.13). The GPS 
displacements are adjusted according the the estimated north-south and east- 
west offsets from Model 2. 

In addition to the GPS data, a priori surface-slip information is added to 
the solution. Surface slip along the Superstition Hills fault (Figure 4.11) 
( Williams and Magistrals, 1989) extends (nearly) the entire length of the 
modeled fault plane. The surface rupture has been incremented into 2.4 km 
segments corresponding to the horizontal dimension of each sub-fault. The 
average slip over each segment is assigned as an a priori slip estimate for the 
surface fault element which it corresponds. Surface rupture along the Elmore 
Ranch fault is confined to the southwestern segment (Figure 4.11). Recall for 
this event that the mapped surface breaks occurred along several nearly 
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Figure 4.10: The residuals (observed minus calculated) for the best-fit 
solution to Model 2. The large residual components in the east-west 
direction are suggestive of measurement error. The residual at MOUN is 
likely attributed to the reset benchmark between surveys. The unusually 
large discrepancy at L589 (and KANE) suggest additional seismic 
deformation not accounted for by the simple uniform slip 
parameterization considered for Model 2. 
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Figure 4.11: Observed surface slip along the Superstition Hills and 
Elmore Ranch faults (dashed lines). The Superstition Hills offsets were 
measured on January 25 and 26, 1988 about 1 month before the GPS 
observations. Decaying afterslip is recorded up to nearly 1 year after the 
earthquake sequence. The Elmore Ranch measurements are the 
cumulative slip from multiple surface breaks across a 10 km wide zone, 
with no recorded postseismic offset after the earthquake. The discrete 
approximation to the surface slip vised to constrain the uppermost sub- 
fault elements in Models 3a and 3b is shown by the solid lines. 
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parallel strands (Figure 4.2). We take the cumulative surface offset for all 
strands [Hudnxit et al., 1989a) averaged over 2.6 km segments along the fault, 
and apply this as a priori slip information for the surface sub-fault elements. 
Where no rupture is mapped (to the northeast), the surface intersecting fault 
partitions are assigned 0 slip. The a priori uncertainty for each surface-slip 
estimate is assumed to be 10 cm. 

The number of independent parameters estimated through singular value 
decomposition depends on the number of singular values k utilized in 
Equation (4.12). A trade-off exists between solution variance and resolution 
[e.g., Menke, 1984). While large k produces highly resolved models, this is at 
the expense of increasing solution uncertainty. Correspondingly, small k 
yields low variance solutions but does not provide detailed resolution. A total 
of 100 sub-fault elements are incorporated into Model 3 (50 for each fault). If 
k — 100 then slip along each sub-fault will be determined uniquely. Because of 
limited geodetic coverage, however, it is practical to consider only the first few 
eigenvectors of the parameter space defined by the geodetic observations. 
Therefore, each solved parameter is a function of some average slip over 
multiple sub-fault elements. This is fundamental property of singular value 
decomposition when used to solve under-determined or poorly-determined 
problems [e.g., Jackson , 1972]. It is necessary to determine the k which 
maximizes the resolution without allowing the solution to become too 
oscillatory or unstable. 

The geodetic moment, solution instability, and model residual calculated 
for different values of k are shown in Figure 4.12. The moment is a function 
of the average slip along the fault planes, while solution instability is 
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determined from the standard deviation of slip for each sub-fault element. 
An instability of 0 (stable) indicates uniform slip along the fault planes, while 
high values indicate an oscillatory or unstable solution. The RMS indicates 
the agreement between model and observation and is calculated by 
RMS «S(°i -c i5/ a i 2 where o, is the observed, c, is the calculated, and a- is 
the uncertainty assigned to the ith observation. 

With surface slip incorporated into the solution (Figure 4.12a), the first 20 
singular values are well constrained by the measured offset along the fault, 
and are influenced little by the geodetic observations. The surface 
measurements reflect the dislocation on the uppermost fault elements, with 
little depth resolution. This is illustrated by the solution for k = 20. The 
calculated moments are significantly less than that for Model 2, presumably 
because the surface offset is not representative of the larger displacement 
along the rest of the fault plane(s). Consequently, in order to estimate the 
slip distribution with depth it is necessary to consider solutions where k > 20. 
After k=»30 the solution becomes very oscillatory as is indicated by the 
increasing instability value. The RMS is significantly reduced beyond k=20 
but only improves marginally with increasing k. The solution fit for k > 20 is 
slightly better than that for Model 2. 

Also evident in Figure 4.12a are the large moment estimates for the 
Elmore Ranch fault, which are almost equal to the computed moments for the 
Superstition Hills fault. This is unexpected considering the latter event 
yielded a significantly greater surface wave magnitude, as well as a larger 
moment estimated from Model 2. It is difficult to distinguish slip between the 
two faults using SVD as is illustrated by Figure 4.13. The displacement 
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Figure 4.12a: The geodetic moment, standard deviation of sub-fault slip 
(instability), and solution RMS calculated for different singular values (k). 
Shown here are solutions constrained by surface slip measurements 
(20 <k <30). 


Constrained 
Surface Slip 


SH 

ER 

Total 



• 228 - 


Figure 4.12b: The geodetic moment, standard deviation of sub-fault slip 
(instability), and sol f S calculated for different singular values (k). 

Shown here are solutic onstrained by surface slip measurements 

(1 <k<10). 
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Figure 4.13: The horizontal slip distribution calculated independently 
for the Superstition Hills and Elmore Ranch faults based on Model 2. 
The shaded region indicates where the horizontal deformation is greater 
than 2 cm for the Superstition Hills fault and 0.5 cm for the Elmore 
Ranch fault. The scale for the Elmore Ranch event is altered to account 
for the smaller dislocation. The deformation pattern is almost identical 
between the two faults, although the displacement magnitudes are larger 
for the Superstition Hills event. This illustrates the difficulty in using the 
GPS measurements to resolve slip between the two faults. 
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pattem for a 130 cm rupture on the northwest trending right-lateral 
Superstition Hills fault is compared with a 30 cm rupture on the northeast 
trending left-lateral Elmore Ranch fault. Although the magnitude 'is different, 
the deformation pattern between the two dislocations is almost identical. As 
a result, the larger slip estimate along the Superstition Hills fault is being 
mapped onto the Elmore Ranch fault plane producing the higher moment. 
This is all the more true considering only the first few model-space 
eigenvectors are utilized, as is necessary since Model 3 is very underdetermined 
(more model parameters than data). Therefore, only linear combinations of 
model parameters are uniquely defined. The similar moments indicated by 
Figure 4.12a suggest slip between the two faults is strongly correlated in the 
solution. 

The estimated seismic slip distribution along the Superstition Hills and 
Elmore Ranch faults for k=23 and k=27 are shown in Figure 4.14a. We 
refer to these solutions as Models 3a and 3b, respectively (Table 4.6). 
Although both faults are partitioned into 50 elements, the fault-rupture is not 
as resolved as the contours suggest since only the first few (non-surface) 
model-space eigenvectors are independently solved. For k = 23 the solution 
suggests fairly uniform rupture along both fault planes. The dislocation may 
be slightly concentrated to the southwest along the Elmore Ranch fault. The 
apparent *bullseye pattern is due to the smoothness constraints requiring the 
slip to tend towards 0 along the lateral edges and lower boundaries (the upper 
boundaries are constrained by the surface slip information). There is little 
difference in the slip distribution for k—21 through k = 26. However, there is 
a noticeable change in the dislocation pattern starting with k =27. While slip 
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along the Elmore Ranch fault still appears fairly uniform, displacement along 
the Superstition Hills fault is concentrated to the northwest and to the 
southeast. This change is significant and is caused by the GPS displacement 
at one station. Recall the large residual for L589 in Model 2 (Figure 4.10). 
This discrepancy is nearly eliminated beginning with k«*27. Therefore, in 
order to satisfy the observed displacement at L589, it is necessary to 
concentrate rupture at each end of the Superstition Hills fault. Of course this 
analysis assumes the observed GPS displacement at L589 is seismically 
generated, and not contaminated by unusually large measurement error. The 
dislocation null near the center of the fault roughly corresponds to the drop in 
fault offset measured at the surface (Figure 4.11). 

Independent solutions are made without constraining the upper sub-fault 
elements by measurements of surface offset (Figure 4.12b). The unconstrained 
moments are generally smaller than when surface slip is incorporated into the 
model. This is because the surface measurements are less than the average 
slip estimate along the fault plane. Without surface constraint, the geodetic 
data are satisfied to a greater degree by slip near the surface; otherwise, it is 
necessary to compensate the small shallow offsets by increased slip at greater 
depths. The seismic slip distribution estimated without surface constraint 
along the Superstition Hills and Elmore Ranch faults for k = 3 and k=7 are 
shown in Figure 4.14b. We refer to these solutions as Models 3c and 3d, 
respectively (Table 4.6). For these unconstrained solutions the slip 
distribution is confined to shallower depths. Slip concentration at each end of 
the Superstition Hills fault is suggested for k » 7, although this division is not 
as pronounced as in Model 3b. We conclude that incorporating measurements 
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Figure 4.14a: Slip distribution along the Superstition Hills and Elmore 
Ranch faults computed from Singular Value Decomposition. Each fault is 
partitioned into the 50 sub-elements indicated by the grid spacing. 
Shown here are solutions for k=23 and k=27 constrained by 
measurements of surface offset. 
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Figure 4.14b: Slip distribution along the Superstition Hills and Elmore 
Ranch faults computed from Singular Value Decomposition. Each fault is 
partitioned into the 50 sub-elements indicated by the grid spacing. 
Shown here are solutions for k = 3 and k = 7 unconstrained by surface slip 
measurements. 
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of surface offset tends to change displacement magnitude by requiring slip at 
greater depths, however, it does not significantly alter estimates of slip 
distribution. 


4.6 Geophysical Implications 
Conjugate Faulting 

The most prominent feature of the Superstition Hills earthquake sequence 
is the conjugate relationship exhibited by near-simultaneous ruptures along 
right-lateral northwest and left-lateral northeast trending faults. In the 
context of the Imperial Valley, the northeast trending structures are termed 
cross-faults [e.g., Hudnut et ai, 1989a]. Conjugate and cross-fault seismicity 
seems to be a fairly typical phenomenon for this region (Figure 4.15), and 
may dictate the strain transfer mechanism between faults. The 1981 
Westmorland earthquake (M L 4.1) is a prime example of cross-fault tectonics. 
The mainshock and aftershock sequence is clearly mapped onto a northeast 
trending lineament. Other examples are associated with the Imperial fault. 
The largest aftershock (M L 5.8) following the 1979 Imperial Valley earthquake 
(M l 6.6) was located near the town of Brawley [. Johnson and Hutton, 1982], 
The focal mechanism and following seismicity suggested left-lateral slip along 
a vertical northeast trending fault. Reilinger and Larsen (1986] found that 
rupture along an identical conjugate structure successfully modeled geodetic 
observations within the Brawley Seismic Zone. A large (M L 5.5) afterahock 
was also recorded near Brawley following the 1940 earthquake [Veumann, 
1942]. Due to the sparsity of seismic data, neither the mechanism nor 
location are precisely determined, although we speculate this event occurred 
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along the same northeast trending feature as the large 1979 aftershock. Of 
historical interest are Imperial Valley earthquake pairs during 1915 (Ml 6.3, 
M l 6.3) and 1927 (M L 5.8, M L 5.5) [BtaL, 1915; Toppozada et ai, 1978]. In 
each case the 2nd shock followed the first by about 1 hour, contrasting with 
the 12 hour interval between the 1987 events. It is not known which fault(s) 
ruptured during these earthquake sequences, but conjugate fault interaction is 
highly probable. 

Rupture on the Superstition Hills fault was almost certainly triggered by 
..^e Elmore Ranch event (occurring 12 hours earlier) suggesting some 
mechanism of stress transfer between the two faults. Figure 4.16 shows the 
normal (<r n ) and strike-shear (cr s ) stress components instantaneously applied 
to the Superstition Hills fault due to a 30 cm left-lateral Elmore Ranch 
dislocation (Model 2). Tension and right-lateral shear are considered positive, 
both tending to induce failure on the rupture plane. Also shown is the 
Coulomb failure stress (cr c ), here given by a c — <x s + fw a , where /i—0.75. 
Positive values indicate stress-loading leading toward shear failure. 

The stress regime necessary for left-lateral rupture along a northeast 
trending fault is identical to that required for right-lateral failure along a 
northwest trending fault. Hence, we can assume that the Superstition Hills 
rupture plane was at or near failure at the time of the Elmore Ranch event. 
The initial shock generated an increase in the Coulomb failure potential along 
the Superstition Hills fault (Figure 4.16), possibly advancing it past its failure 
threshold. This is seen mostly as a combination of reduced compressive 
normal stress (earthquake inducing) countered by left-lateral shear 
(earthquake inhibiting). The increase is maximized along the northwest 
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Figure 4.15: Known and/or potential conjugate/cross-fault seismic 
episodes in the Imperial Valley since 1900. Seismic release on (left-lateral) 
northeast trending structures was observed in 1979, 1981, and 1987. 
Earthquake pairs or mainshock/aftershock sequences suggestive of 
conjugate faulting were observed in 1915, 1927, 1940. This suggests 
conjugate/cross-fault interaction is typical for the Imperial Valley. 









normal stress 
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Figure 4.10: Strike-shear (right-lateral positive) and 
change (dilatation positive) induced on the Superstition Hills fault due to 
a 30 cm left-lateral Elmore Ranch dislocation (Model 2). Also shown is 
the Coulomb failure stress change, where positive values indicate an 
increased potential for rupture (earthquake inducing stress). The 
northwest third of the 1987 Superstition Hills rupture plane underwent a 
stress change tending it toward failure with the maximum change 
calculated in the epicentral region near the intersection with the Elmore 
Ranch fault. There was no rupture to the northwest where the Coulomb 
failure stress was negative (reduced earthquake potential). The 

magnitude of the stress change (in bars) is in the range of typical 
earthquake stress drops. 



Depth (km) Depth (km) Depth (km) 


- 243 - 


0 
5 
10 
0 
5 

10 
0 
5 
10 

-20 0 20 
Distance Along Strike (km) 






1-2 


5-10 


> 10 









- 244 - 


boundary of the rupture plane, near the nucleation point of the second event. 
Presumably rupture began where the applied stress was greatest and then 
propagated to the southeast. Northwestward rupture is prohibited because 
the increase in compressive forces tends to inhibit shear failure along this 
segment of the fault. The magnitude of the Coulomb stress increase near the 
Superstition Hills epicentral zone (~ 10 bars) is comparable to typical 
earthquake stress drops. 

The one to several hour delay recorded between events during observed 
and suspected conjugate episodes in the Imperial Valley is significant from an 
earthquake failure perspective. Shown in Figure 4.17 are potential scenario’s 
for earthquake ruptures involving conjugate-mainshock interaction, such as 
that observed for the Superstition Hills events. We assume faults fail by an 
undefined mechanism when they are at or above some critical str ess level. 
The regional strain acting over several years brings a given fault near this 
critical failure point. A stress increase is induced along part of the fault plane 
due to rupture on a conjugate structure (e.g., Figure 4.16), which may or may 
not be sufficient to push the stress state past its critical threshold. In the case 
of Earthquake 1 (Figure 4.17a), the stress change caused by the conjugate 
event is not enough to induce failure. Some form of time dependent stress 
transfer onto the fault is active and eventually the critical level is reached. 
Such a mechanism involving postseismic viscous creep along the Elmore 
Ranch fault has been suggested for the 1987 Superstition Hills sequence [Given 
and Stuart , 1988). If this scenario is valid we would equally expect failure 
modes such as that indicated by Earthquake 2. Here the instantaneous stress 
applied to the fault from the conjugate event pushes the stress state past the 
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critical level and rupture is immediate. In this case failure along the two 
perpendicular fault planes will occur simultaneously. However, this behavior 
is not observed in the Imperial Valley. Conjugate episodes characteristically 
have been separated by one to several hours. This suggests that the critical 
stress level can be exceeded without immediate failure. Therefore, some time 
dependent mechanism must be active on the fault plane, as opposed to 
additional stress transfer through the crust. We loosely refer to this as stress 
corrosion" (Figure 4.17b) [e.g., Das and Scholz, 1981]. In the case of 
Earthquakes 3 and 4, it is suggested that the critical stress level must be 
exceeded for a period of one to several hours before failure occurs. Hxidnut et 
al. (1989b] proposed fluid diffusion as an alternate mechanism, whereby the 
effective normal stress was reduced (made more positive) due to pore-fluid 
infiltration into the rupture plane, thus increasing the Coulomb failure stress. 
This process still involves action on the fault plane rather than stress transfer. 
Regardless of cause, the temporal and geometric relationship exhibited by the 
conjugate fault interaction is seemingly typical of Imperial Valley tectonics 
and is thus an important factor for the potential prediction of large 
earthquakes and aftershocks. 

Moment and slip distribution 

The geodetic (GPS) source parameters for the Superstition Hills and 
Elmore Ranch earthquakes are listed in Table 4.6 and Figure 4.14. The 
seismic moment is best defined by Model 2, while the slip distribution is best 
expressed by Models 3a and 3b. Model 1 is constrained with minimal station 
coverage and Models 3c and 3d do not include the additional information 
supplied from surface slip measurements. The GPS observations are directly 
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Figure 4.17: Schematic of potential earthquake failure processes in the 
Imperial Valley, a) Earthquake failure occurs after some critical stress is 
reached, b) Earthquake failure occurs following a time dependent delay 
after critical stress is exceeded. 
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proportional to the combined effect of the Elmore Ranch and Superstition 
Hills events, although we have attempted to resolve slip between each fault 
plane. The calculated parameters are a function of the coseismic oflset, as 
well as 3-4 months of postseismic slip (plus 1.5 years of preseismic movement, 
if any). The average Elmore Ranch dislocation is about 30 cm (left-lateral) 
with fairly uniform distribution along the fault plane. In the case of the 
Superstition Hills fault, the average slip is estimated at 130 cm with 
concentrated deformation along the northwest and southeast sections of the 
fault. Because the GPS sampling frequency is so low (years), the calculated 
source parameters should contain the total coseismic moment release, which 
includes several months of postseismic slip. 

In Table 4.7 the GPS moments are compared with estimates made 
through seismic and other geodetic studies. Forward and inverse models 
using teleseismic [Dziewonski et a/., 1989; Bent et ai, 1989; Sipkin, 1989; 
Hwang et al. , 1990) and strong-motion recordings [Frankcl and Wennerberg , 
1989; Wald et ai, 1990] are used to constrain source parameters, as well as 
investigate complexities of the Superstition Hills rupture process. The 
teleseismic moments agree fairly well with the GPS estimates, while the strong 
ground motion data yield significantly lower results. The high-frequency 
strong-motion measurements are dominated by energy around 1 second and 
conceivably miss a sizable portion of the long-period energy release recorded 
with the teleseismic and GPS data. Hence, the near field seismic solutions 
may underestimate the total moment release. 

Geodetic measurements from Pinyon Flat observatory are used to 
constrain planer and curved dislocation models for the Superstition Hills and 
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Table 4.7 Moment Comparison 


Method 

Moment (xlO 26 dyne-cm) 
SH ER Total Ratio 

Reference 

GPS (Model 2) 

9.4 

2.3 

11.7 

4.1 

This Study 

Teleseismic 

7.2 

1.4 

8.6 

5.1 

Dziewonski et ai [1988] 

Teleseismic 

10. 

2.3 

12. 

4.3 

Sipktn [1989] 

Teleseismic 

10.8 

2.7 

13.5 

4.0 

Bent et ai [1989] 

Teleseismic 

8. 




Hwang et al. [1990] 

Strong Motion 

5.2 




Wald et ai [1989] 

Strong Motion 

1.8 




Frankell and Wennerberg [1989] 

Pinyon Flat (Planar- A) 

3.7 

0.8 

4.3 

4.6 

Agnew and Wyatt [1989] 

EDM 

9.3 




LUomki and Savage [1988] 


SH - Superstition Hills fault 
ER • Elmore Ranch fault 
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Eimore Ranch faults \Agnew and Wyatt, 1989). The data are obtained from 
long-base strain and tilt-meters, as well as a borehole dilatometer. The best- 
fit planer models to all observations (Table 4.7) are significantly lower than 
those calculated with the GPS and teleseismic data, although a 70 % moment 
increase for the Superstition Hills fault is obtained when the strainmeter data 
are excluded. The low moment estimate may be due to a number of factors 
[Agneru and Wyatt, 1989): 1) measurement quality, particular with the 
strainmeter, 2) rheologic differences between Superstition Hills and Pinyon 
Flat, and 3) strainmeter-dilatometer sensitivity to the nodal deformation 
plane roughly on azimuth with the observatory. 

Geodolite observations of the Salton Trough EDM network were made in 
early December (1987), several days after the two large earthquakes [Lisowski 
and Savage, 1988); the last previous occupation was in January, 1987. Simple 
dislocation models with 40 cm left-lateral slip along the Elmore Ranch fault 
and 120 cm right-lateral slip along the Superstition Hills fault best-fit the 
observations. The estimated moment for the Superstition Hills event (Table 
4.7) is comparable to that obtained with the GPS displacements. 

The discrepancies in Table 4.7 are largely attributed to the alternate 
methodologies, observations, and parameters used to constrain each model. 
However, for those calculations which include moment estimates for both the 
Superstition Hills and Elmore Ranch events, the ratio between the two 
ruptures is fairly constant. This illustrates an internal consistency with each 
method. More importantly, it suggests that postseismic slip along the 
Superstition Hills fault is probably confined to the shallow segment of the 
rupture plane. Since seismic activity on the Elmore Ranch fault essentially 


- 251 - 

ceased after the 2nd main event, if postseismic slip were occurring in mass 
along a large fraction of the Superstition Hills rupture plane, the GPS 
moment ratio would be considerably larger. 

While the epicenter and aftershock sequence for the Superstition Hills 
event were concentrated along the northwestern portion of the fault, strong 
ground motion, teleseismic, and surface offset data suggest significant moment 
release on the southern section of the Superstition Hills fault [Wald et ai, 
1989; Bent et ai , 1989; Hwang et ai , 1990; Williams and Magistrate , 1989). 
An exception is the strong ground motion study of Frankel and Wennerberg 
[1989] where slip is confined to the northwest. However, their low 
Superstition Hills moment (Table 4.7) suggests this analysis may be strongly 
susceptible to the high-frequency content of the data, indicating rupture along 
the southeast segment was dictated primarily by low-frequency energy release. 
The GPS data also reveal dislocation along the southeastern segment of the 
fault, and further suggests a displacement null near the faults mid-section. 
This slip deficiency may be related to the surface offset drop observed along 
the center of the fault (Figure 4.10). 

Deformation across the Imperial Valley 

The 1986-1988 GPS station displacements indicate a significant 
component of deformation across the Imperial Valley not associated with the 
1987 Superstition Hills earthquake sequence (Figures 4.4 and 4.8). This 
motion is attributed to plateboundary deformation due to the relative 
velocity between the Pacific and North American plates. From empirical 
evidence provided by Salton Trough EDM observations between 1972 and 
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1087, these non-seismic movements were modeled as 3.45 cm/yr differential 
velocity across the valley. However, after removing the seismic deformation 
predicted with a preliminary model (Model 1, Table 4.6) the calculated GPS 
displacements average 6.0 cm/yr (Figure 4.8), significantly larger than that 
obtained with the trilateration measurements over the last two decades. This 
GPS velocity is also larger than the 4.7 cm/yr plate motion predicted from 
global models (for Imperial Valley coordinates) (. DcUct a ef «/., 19g0 |. Most 
likely the GPS rate is an overestimate, especially considering it contains large 
measurement uncertainty due to poor 1086 data quality, and is heavily 
influenced by seismic effects from the 1087 earthquake sequence. It is possible, 
however, that the 1086-1088 GPS rat. could represent accelerated’ 
deformation. In fact, GPS observations between 1088 and 1080 indicate 5.2 ± 
0.0 cm/yr displacement across the Imperial Valley [iareen and RaMnger, 
1000]. Accelerated deformation is not without precedence. Triangulate 
observations suggest a rate of 6.2 cm/yr between 1041 and 1054, although 
this is attributed to postseismic deformation following the 1040 Imperial 
Valley earthquake. Increased deformation following the 1070 earthquake is 
not observed in the EDM observations [Savage et al„ 1086|. Additional GPS 

measurements are necessary in the Imperial Valley in order to ascertain the 
current deformation rate. 


4.7 Conclusions 

Station movements computed from 4 Imperial Valley GPS campaigns 
indicate large crustal displacements during the periods 1986-1988 and 1986- 
1990. Much of the deformation is attributed to the 1987 Superetition Hills 
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earthquake sequence. Ten sites near the seismic rupture zone are displaced at 
least 10 cm, although the GPS observations contain large uncertainties due to 
poor data quality from the initial (1986) survey. This is the first occurrence 
of a large earthquake within a preexisting GPS network. 

The 1987 earthquake sequence is initially modeled as uniform offsets along 
rectangular dislocations in an elastic half-space. The best-fit model to the 
GPS observations requires 130 cm right-lateral slip along the northwest 
trending Superstition Hills fault and 30 cm left-lateral motion along the 
conjugate northeast trending Elmore Ranch fault. The slip distribution along 
each fault is investigated by partitioning the rupture planes into 50 sub- 
elements and utilizing Singular Value Decomposition to estimate the slip 
along each sub-fault. Measurements of surface offset are used to constrain the 
shallow elements of the fault plane. The estimated slip distribution along the 
Elmore Ranch fault is fairly uniform. Slip along the Superstition Hills fault 
appears to be concentrated to the northwest and the southeast with a 
displacement drop indicated near the faults midsection. There is some 
evidence that postseismic slip along the Superstition Hills fault was 
concentrated near the surface. The estimated moments are 9.4 x 10 2S dyne- 
cm and 2.3 x 10 25 dyne-cm for the Superstition Hills and Elmore Ranch 
faults, respectively, which are consistent with moments obtained from 
teleseismic data. 

In addition, the 1986-1988 GPS observations suggest non-seismic 
movements across the Imperial Valley of up to 5.9 cm/yr. These secular 
displacements are attributed to pi ate- boundary deformation due to the 
relative motion between the North American and Pacific plates. The observed 
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rate is probably an overestimate, however, as it is heavily influenced by 
unmodeled seismic effects and measurement error. Regardless, the observed 
seismic and secular deformations clearly emphasize the importance of future 
GPS study in the Imperial Valley. 
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ABSTRACT 


Resolution of both the extent and mechanism of lateral heterogeneity in the upper 
mantle constrains the nature and scales of mantle convection. Oce^egrons I “* f of 
particular interest as they are likely to provide our closest glimpse at the patterns of 
temperature anomalies and convective flow in the upper mantle because of fhe^ yowij age 
Stole crustal structure relative to continental regions. Our objectives in this thesis are 
to determine lateral variations in the seismic velocity and attenuation [Structure of the 
lithosphere and asthenosphere beneath the oceans, and to combine these seismotogical 
observations with the date and theory of geoid and bathymetry anomalies in order to test 
an^ improve current models for seafloor spreading andmantl^onvecuon^ ^corteenffate 
on determining variations in mantle properties on a scale of about 1000 km, comparable to 
the thickness of the upper mantle. Seismic velocity, geoid, and bathymetry anomalies are 
all sensitive to variations in upper mantle density, and we formulate inversions to combine 
quantitatively these different date and search for a common origin. Vanauot^m mantle 
density can be either of thermal or compositional on gin and are presumably related 
mantle convection and differentiation. 

Bv means of a large date base of digital seismograms and waveform cross-cwrelation 
and sStio techniques, we have measured SS-S differential travel time residuals and 
differential attenuation in order to determine lateral variations in upper mantle str ^“]' 
beneath the Mid-Atlantic Ridge and East Pacific Rise. Differential travel times of such 
chases as SS and S with identical source and receiver have the advantage that residuaisare 
UlSly to be dominated by contributions from the upper mantle next 
of the reflected phase (SS). Under this assumption, differential SS-S travel time residuals 
are mapped at the SS bounce points as a means of delineating lateral variations m mande 
structure^After removing the signature of lithosphere age, we find evidence for long- 
wavelength variations in SS-S residuals along the Mid-Atlantic Ridge. The dormnan 
wavelength of these variations is 1000 to 2000 km. These travel time anomalies correlate 
qualitatively with along-axis variations in bathymetry and geoid height We formulate a 
joim inversion of travel time residual, geoid height, and bathy^ under the assumption 
that all arise from variations in upper mantle temperature or buUe compoanon 
(parameterized in terms of Mg#). The inversion employs geoid and topography kernels 
which depend on the mantle viscosity structure. Inversion for temperature perttebations 
2S ^SsTood fits to travel time and geoid data. The fit to topography,which is 
jjfcely dominated by unmodeled crustal thickness variations, is not as good. The in ^®? lonS 
for temperature favor the presence of a thin low viscosity layer in the upper mantle and 
temperature perturbations concentrated at depths less than 300 ^- Com ^ s m° n 
variations alone are unable to match the travel time and geoid or bathymetry date 
simultaneously. A joint inversion for temperature and composition provides ; good fits to 
both geoid and travel time anomalies. Temperature variations are ± 50 K and 
compositional variations are ± 0.5-3 % Mg# for models with the temperature vmauons 
uniformly distributed over the uppermost 300 km and the ^mposioon^vmanons eitiier 
distributed uniformly over the same interval or concentrated at shallower depths. The 
magnitudes of these variations are consistent with the chemistry and geothermometry of 
dredged peridotites along the Mid-Atlantic Ridge. 

nifWntial travel times of SS-S pairs in the east central Pacific show several 
differences from the north Atlantic. The most obvious difference is that the travel time 
residuals are significantly larger than in the Atlantic, even at a fixed age. The travel tune- 
age relation is weaker in the Pacific, although this may be pamaily atmbuteble to *e fact 
SSt we have not sampled a large range of plate ages in the eastern Pacific. In the Atlantic 
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our results are not consistent with the presence of a simple pattern of tmtsotropy^ 

while in the Pacific the data are consistent with the presence of weak anisotropy in the uppe 
mantle It has been suggested that anisotropy may be more pronounced at fast spreading 
^sto a^spSgrares both in ft? lithosphere (due to a rate dependence of the 
mechanism for orienting olivine crystals in the lithosphere) and the asthenosphere (because 
the asthenospheric flow 8 beneath fast moving plates is likely to lake the form of a progressive 
simple shear which can produce a lattice preferred orientation of olivine crystals), and our 
Site SSS this suggestion There is substantial am^mtyin ouranisom>py 
measurements for the Pacific, however, due to a poor sampling of azimuthsjo mat it is 

also possible that lateral heterogeneity rather than a ^ uth ‘f^ 1 ^ 
observed azimuthal pattern. Sampling at a more uniform distribution of ^nmthsj should 
make this result less ambiguous, and as more seismic stations are deployed at new 
geographic locations our chances of resolving mis issue will improve. 

Inversion of travel time residuals, geoid, and bathymetry data for jhe < ®£ t ^J a PaClflC 
indicates that compositional variations alone are inadequate to match all of the data 
rimS^sly^Sar to our results for the north Atlantic. Tempm^ van^ons^, 
howevS, produce significant variance reduction. The inversion sduncm indic^eexcess 
temperature in me vicinity of me Galapagos hotspot in me range 50 - 150 K. Furthe , 
anal>Ss1sneeded to determine me effects of subduction zone structure and possible crustal 

thickening in me eastern Cocos plate region. 

As a complement to me study of travel times, we have measured SS-S differential 
attenuation in me north Atlantic region. Mapping seismic Q in me upper mantle is an 
important tool for assessing mechanisms of lateral heterogeneity because tiw a^nuanonof 
seismic waves is sensitive to variations in temperature and to parttal melon g. Differe 
SSi^ Svely correlated with SS-S travel time residual. Bom differential 
attenuation and travel time residual decrease with increasing seafloor age. The age 
dependence of SS-S travel time residual can be explained entirely by me cooling of the 
oceanic lithosphere i e., contributions from me asthenosphere or from a mantle melt 
faction are not required. On me assumption mat plate cooling also doimmtes the vmauon 
of differential attenuation with age, we derive an empirical Q* -temperature relation for t 
"SfespE Se variation of Q > with temperature that we derive is not as strongly 
dependent on temperature as mat observed in laboratory studies. Systematic long 
wsn^leneth ( 1000-6000 km) variations in upper mantle differential attenuation are ev ^d 
alone theaxis of me Mid- Atlantic Ridge. These variations correlate approximately with 
feng- wavefen^h Variations i nshear wfve travel time residuals and are attributed to along- 

axis differences in upper mantle temperature. 
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ABSTRACT 

Utilizing a digital data base of over 500 seismograms and a waveform cross-correlation 
technique, we have measured SS-S differential travel time residuals as a means to examine the 
shear wave velocity structure in the vicinity of the Mid-Atlantic Ridge. Differential travel times 
of such phases as SS and S with identical source and receiver have the advantage that residuals 
are likely to be dominated by contributions torn the upper mantle near the surface bounce point 
of the reflected phase (SS). Under this assumption, differential SS-S travel nme residuals are 
mapped at the SS bounce points as a means of delineating lateral variations in mantle structure. 
After removing the signature of lithosphere age, we find evidence for long-wavelength variations 
in SS-S residuals along the ridge. The dominant wavelength of these variations is 1000 to 2000 
km. These travel time anomalies correlate qualitatively with along-axis variations in bathymetry 
and geoid height We formulate a joint inversion of travel time residual, geoid height and 
bathymetry under the assumption that all arise from variations in upper mantle temperature or 
bulk composition (parameterized in terns of Mg#). Hie inversion employs geoid and 
topography kernels which depend upon the mande viscosity structure. Inversion for thermal 
perturbations alone provides good fits to travel time and geoid data. The fit to topography, 
which is likely dominated by unmodeled crustal thickness variations, is not as good. The 
inversions for temperature favor die presence of a thin low viscosity layer in the upper mantle 
and temperature perturbations consented at depths less than 300 km Compositional variations 
alone are unable to match the travel time and geoid or bathymetry data simultaneously. A joint 
inversion for temperature and composition provides good fits to both geoid and travel time 
anomalies. Temperature variations are ± 50 K and compositional variations are ± 0.5-3 % Mg# 
for models with the temperature variations uniformly distributed over the uppermost 300 km and 
the compositional variations either distributed uniformly over the same interval or concentrated a. 
shallower depths. The magnitudes of these variations are consistent with the chemistry and 
gcothermometry of dredged peridotites along the Mid-Atlantic Ridge. 
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INTRODUCTION 


Seismic velocity and density of upper mantle material are expected to be functions of 
temperature and composition. The delineation of long wavelength variations in these physical 
properties thus provide important constraints on mantle convection, crust-mantle differentiation, and 
mantle chemical heterogeneity. In this study we determine lateral variations in upper mantle 
temperature and composition along the Mid-Atlantic Ridge through the combined inversion of shear 
wave differential travel times, geoid height, and bathymetric depth anomalies. 

The advent of seismic tomography has led to a number of three-dimensional maps of lateral 
variations in seismic velocity in the upper mantle, and several such models of the north Atlantic 
region have been developed, both as parts of global studies [e.g., Woodhouse and Dziewonski , 
1984; Nakanishi and Anderson, 1984; Tanimoto, 1990] and through regional investigations of 
long-period surface waves [e.g., Honda and Tanimoto , 1987; Mocquet et al, 1989; Mocquet and 
Romanowicz, 1990]. With surface wave methods each wave samples the average vertical variation 
in upper mantle structure along its path, but because of the long wavelengths involved the inversion 
of phase or group velocity from many paths tends to smooth out lateral variations. Body wave 
travel times can provide independent information about upper mantle heterogeneity at potenually 
shorter horizontal scales than surface waves can resolve, and progress has been made in the 
determination of lateral heterogeneity in the North Atlantic through the use of both differential and 
absolute travel times of body waves [Kuo et al., 1987; Grand, 1987, 1989]. 

The travel times used in this study are differential times of the body wave phase pair SS-S. 
Differential travel times of shear wave pairs are well suited to the study of upper mantle 
heterogeneity [Sipkin and Jordan, 1976, 1980; Stark and Forsyth, 1983; Butler, 1979; Kuo et al, 
1987; Woodward and Masters, 1991] and have the advantage that source and receiver effects are 
approximately common to both phases and are thus largely eliminated by differencing. Under the 
assumption that the lower mantle is relatively homogeneous and that the portions of the ray paths in 
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the upper mantle are steep, the differential travel time anomaly is associated with upper mantle 
heterogeneity in a small volume centered beneath the surface bounce point of the reflected (SS) 
phase. This technique is thus well suited to the investigation of horizontal variation, m structure, 

but the resolution of variations with depth is poor. 

Oceanic bathymetry and geoid height data are sensitive to variations in mantle density at 

depth. Such variations can be either thermal or compositional in origin and, like seismic velocity, 
are presumably related to mantle convection and differentiation. Geoid (or gravity) and topography 
have become the most commonly used tools for mapping out and constraining models of upper 
mantle convection [e.g., Anderson et al., 1973, McKenzie and Bowin , 1976; McKenzie, 1977; 
McKenzie et al., 1980, Parsons and Daly , 1983; Buck and Parmentier , 1986; Craig and McKenzie, 
1986]. In addition, measurement of the admittance (the spectral ratio of geoid to topography) has 
been widely utilized to estimate the depth and mode of compensation of oceanic swells and plateaus 
[e.g.. Watts et al. , 1985, Cazenave et al. , 1988; Sandwell and MacKenzie, 1989; Sheehan and 
McNutt, 1989]. Several workers [Dziewonski et al., 1977; Nakanishi and Anderson, 1984; 
Tanimoto and Anderson, 1984; Stark and Forsyth, 1983; Dziewonski, 1984; Kuo et al., 1987] 
have noted coirelations of geoid and travel time (or velocity structure) at a number of different 
wavelengths, although only a few [Richards and Hager, 1984; Hager and Clayton, 1989] have 
combined observational seismology with geoid anomalies in a quantitative and dynamically 

consistent manner. 

In this study we present the first formal inversion of geoid, depth, and travel time anomaly 
data for lateral variations in upper mantle temperature and composition along the Mid-Atlantic 
Ridge. Given a distribution of temperature or density perturbations in the upper mantle, the 
forward problem of calculating differential travel time, geoid, and depth anomalies is 
straightforward. This forward problem forms the basis for a joint linear inversion of these three 
types of observations under the assumption that all arise from parameterized long-wavelength 
variations in upper mantle temperature or composition. Results of a set of inversions carried out 
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under different assumptions regarding the depth extent of lateral heterogeneity and the mande 
viscosity structure are compared with other constraints on variations in mantle temperature and 

degree of melt removal. 


MEASUREMENT OF DIFFERENTIAL TRAVEL TIMES 

The seismic data used in this study consist of long-period S and SS phases obtained from the 
Global Digital Seismic Network (GDSN) [Peterson et al„ 1976; Peterson and Hast, 19821; the 
Network of Autonomously Reconiing Seismographs (NARS). a linear broadband array in western 
Europe [Note, and Vlaar, 19821; and several broadband stations from the global GEOSCOPE 
network [ Romanowicz et at., 1984). A list of stations used in this study is presented in Table 1. 

We use only transversely polarized (SH) seismograms (totaled from N-S and E-W components) to 
avoid interference from die SKS phase and contamination from P-SV conversions at the base of the 
crust and other near-surface discontinuities. Recent work by Gee and Jordan [1989] suggests that 
navel times depend on the frequency band used in the analysis. In order to maintain a self- 
consistent data set for our study, additional processing is applied to data from the NARS and 
GEOSCOPE arrays in older to mimic the instrument response of the longer period GDSN stauons. 
This processing allows us to measure travel times from a set of seismograms that all have 
essentially die same frequency response. Data from the NARS and GEOSCOPE arrays are 
decimated (with a low-pass antialiasing filter) to a common sampling interval of 1 s. The data are 
further filtered using a noncausal 3-point Butterworth filter {Rader and Gold, 1967] with a 
frequency bandpass of 0.01 - 0.20 Hz. This additional filtering gready improves the signal-to- 

noise ratio of the SS phase. 

A waveform cross-cotrelation method is utilized to determine the differential travel time 
between the phases S and SS {Butler, 1979; Stark and Forsyth, 1983; Kuo et al„ 1987], The 
procedure involves the construction of a "synthetic" SS pulse from S and the evaluation of the 
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cross-correlation function between the real and synthetic windowed SS phases (Frgure 1). The 
synthetic SS pulse is created Iron, S in tire following manner. The S pulse is windowed and 
attenuated (with attenuation parameter «- 3 s) iGrandandHelmberger, 1984; Kuo cal, 1987] to 
account for tire additional time SS travels in tire mantle, and then a rt/2 phase shift (Hr, her. 
transform) is applied to tire attenuated S pulse to simulate tire frequency-dependent phase shrft 
which meSS wave undergoes a. an internal caustic { Choy and Richards. 1975]. H,e drfferenn 
dme is obtained from the peah of the cross correlation between tire synthetic SS consumed from 
tire S wave and tire real SS. The residual SS-S times are obtitined by subtracting the obsrrved 
differential time from that predicted by the PREM Ear* mode. and Anderson, ,981, 

and correcting for Ear* eltipticity [Dzieu^ and Ohr «. >976, and SS bounce point 
bathymerry. Our convention is drat negative residuals am indicative of either early SS or latir S. 

Constant window lengths of 120 s are used for both the S and SS phases. In general, the 
observed differential travel times vary by as much as 1 s dependtng on how S and SS are 
windowed. Our modelling with synthetic seismograms indicates that emphasizing tire onset of tire 
SS waveform can lead to bias for bounce points in areas of oceanic sediments. The effect of 
sediments a, long periods is to prmiuce precursory arnvals from reflections a. tire base of the 
sediments and late arrivals from waves which travel through the low-veloctty sediments an 
reflected at the crus.-wamr interface. Tte ne, affecti after convolving tire eras*, response witi, tire 
long-period GDSN instrument response, is that the time center of the SS phase is effecti e y 
unchanged bu, tire pulse is bmadened both a, tire front and a, tire bach. In our procedure the use of 

• • rvnrirsv 9S n »lse should yield differential travel times that are little 
a constant window containing the entire SS puls y 

affected by the presence of sediments. 


DATA 


The north Atlantic is an 


ideal area for conducting a 


differential travel time study in terms of 




7 


ihe geographic distribution of available events and stations at suitable distances. The range in 
source-receiver separation was taken to be 55’ to 86* to ensure separation of S and ScS at the longer 
distances and to avoid triplication in SS at shorter distances. The SS and S phases bottom from 
about 670 km to 2300 km depth. We performed a search over all earihquakes in the Harvard 
centroid moment tensor (CMT) catalog for die years 1977-1987 mewonskietaL 1981; 

Dziewonski and Woodhouse, 1983] and over all GDSN, NARS. and GEOSCOPE digital seismic 
stations in order to find event-station pairs of the proper epicentral distance which provide SS 
bounce points in the North Atlantic region. Epicenters were obtained from the Bulletin of the 
International Seismological Centre (ISC) for events occurring before 1987 and from the 
"Preliminary Determination of Epicenters" of theU.S. National Earihquake Information Service 
(NEIS) for events occurring in 1987. The final distribution of sources and stations used to measure 
SS-S differential travel times is shown in Figure 2. The majority of data in this study comes from 
records of equatorial fracture rone earthquakes at North American and European stations, north and 
central Atlantic events at North American stations. Central American events at European stations, 

and Mediterranean and European earthquakes at North American stations. 

This search yielded over 2000 event-station pairs with the proper epicential separation. After 
winnowing the list because of station inoperation, poor signal to noise ratio for the phases of 
inrerest, and interfering events, tire final data se, consists of nearly 500 SS-S differential nave, time 
residuals with bounce points in the north Atlantic (Figure 3). Uncenainties are determined for each 
measurement following the procedure outlined in Appendix A. 


RESULTS 


We interpret the variations in SS-S differential travel times in terms of lateral veloctty 
variations within tile eras, and upper mantle beneath die surface reflection points of the SS wave 
path. Kuo e, at. [1987] and Woodward and Masters [1991] tested the validity of this assumption 
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by plotting absolute S and SS residuals against SS-S residuals. They found that S and SS-S 
residuals are unconelated while SS and SS-S residuals are strongly correlated, indicating that the 
assumption is justified. The validity of this assumption is further supported by the strong 
correlation of SS-S times with surface tectonic features in the vicinity of the SS bounce point. The 
residuals are further interpreted in terms of such upper mantle processes as lithospheric aging, flow- 
induced anisotropy, and along-axis heterogeneity in mantle structure. 

Lithospheric Aging 

Cooling and thickening of the lithosphere should yield a tendency toward an increase in 
seismic velocity with increasing lithospheric age. A linear regression experiment was performed to 
examine the correlation of the SS-S residuals with seafloor age. A gridded map of seafloor ages 
was constructed for the north Atlantic from the magnetic anomalies of Klitgord and Schoulen 
[1986] and ages assigned according to Kent andGradstein (1986) and Klitgord and Schoulen 
[1986], The isochrons of Sclater el al. [1981] were used in a few regions which were not covered 
by the Klitgord and Schoulen [1986] data set. To obtain a representative age value for the regton 
spanning approximately one horirontal wavelength of the incident (SS) wave, an average seafloor 
age was estimated for a 1' x 1* box centered on each SS bounce point. To reduce scatter, 
measurements whose bounce point depths differed by more than 2500 m from the depth predicted 
by the Parsons and Sclater [1977] plate cooling model were excluded from the final age regression. 
Although each SS wave samples the upper mande at a finite range of lithosphere ages, we expect 
dutt the different travel time anomalies contributed by the SS path segments on the younger and 
older sides of the bounce point approximately cancel so that the age at the SS bounce point is 
appropriate to the associated SS-S residual. 

The SS-S residuals for the north Adamic are consistent with the expectation of an increase in 
seismic velocity with seafloor age. For bounce points between O’ and 60'N latitude, the coefficient 
derived by linear regression of residual with square root of age is -0.68 ± 0.08 s My W from 0 to 
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100 My, with a linear cotrelation coefficient of -0.85 (Figure 4). However, residuals from 60 - 
90-N do not seem to be strongly correlated with lithospheric age. This may be due to the fact that 
this area is more tectonically complicated than 'normal' oceanic lithosphere [e.g., White, 1988; 
Zehnder and Mutter, 1990], includes several ridge jumps, is in close proximity to continental 
regions, and does not closely follow the age-depth relation of Parsons and Sclater [1977], 

Compared with the residuals for OWN, those from 60-90'N are anomalously negadve at young 
ages and anomalously positive at older ages. The slope of SS-S residual vs. square root of age for 
data from 0 to 60'N is smaller than that inferred from S delays of intraplate earthquakes in the 
Atlantic by Duschenes and Solomon 1 1977] (two-way S delay = - 1 .2 s My-'*) and that reponed by 
Kuoetal. [1987] (-1 s My m ). It is larger, however, than the global average obtained by 
Woodward and Masters [1991] (-0.51 s My>«). We find that the residual-age relation is not 
constant over the entire north Atlantic, so that some of these variations in slope may reflect real 
geographic differences. 

We may compare the variation of SS-S residual versus age with that due only to lithospheric 
cooling. For a lithospheric structure given by the plate cooling model of Parsons and Sclater 
[1977], we may convert temperature variations to differences in shear velocity v s by adopting a 
value for 8vs/8T, which we take to be uniform and equal to -0.6 m/s K’ 1 [McNutt and Judge , 

1990]. For a horizontal slowness typical of the teleseismic S and SS waves of this study (0.1375 
s/km), the slope of the line best fitting the SS-S travel time delay versus age given by the plate 
cooling model over 0-100 My is then -0.73 ± 0.07 s My 1/2, a result indistinguishable from the 
observed slope. This agreement indicates that the dependence of travel time residual on plate age 

can be explained entirely by lithospheric cooling. 

The trend of the travel time residual versus lithospheric age relation changes at about 100 My 
After 100 My, the residuals appear to flatten out (Figure 4), in the same sense as the plate cooling 
model of Parsons and Sclater [ 1977], Such a pattern may reflect the unmodeled effect of increased 
sediment or crustal thickness, or, as suggested by Parsons and Sclater [1977] may be partially the 
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rcsu n of second convection which supplies heat to dro base of the plate a, older ages. To avotd 
possible biases associated with any of these effects we shall testnc oar analysis to data wtth bounce 
points on lithosphere iess than 100 My. To Ieoh for other systematic variations » the restdaals, we 
correct for age by removing the linear relation shown by the solid tine in Figure 4. This correction 
is effectively a normalization of residuals to 22-My-old lithosphete (the zero crossmg of the 

regression line). 


Anisotropy 

Another systematic velocity variation that has been suggested as a possible conmbutor to 

residual SS-S travel times is azimuthal anisotiopy. Kuo « * OT P “ 0 " 

detail and concluded that alignment of olivine cystals in the asthenosphete created a stgntftcant 
panem of azimuthal anisouopy in SS-S residuals measured in the Atlantic region. We have also 

searched for evidence of azimuthal anisotropy with our data set. 

Backus [1965] and Crampin [1977] demonstrated, from the general form of body wave 

anisotropy in a weakly anisouopic medium, .ha, the linear form of the azimuthal variation of 
velocity is given by 

V 2 = Aq + A x cos 20 + A 2 sin 20 + A 3 cos 40 + A 4 sin 40 (1) 

where V is the tody wave velocity, tit. A„ are linear functions of the eiastic moduli, and 6 is an 
azimuth, defined for our problem by the angle between the groat cirole path and tiro drrocnon to 
geographic nortit measured at the SS bounce point Equation (1) was furtiter simphfied by Kuo e, 
al. 119871 and parameterized in terms of travel time residuals: 


R = Ro + R 1 cos 20 + R 2 sin 20 + R 3 cos 40 + R 4 sin 40 


( 2 ) 
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where R is the travel time residual and the R„ are constants. By fitting a function of this form to our 
age-corrected measurements we can determine if our data are consistent with the presence of 

anisotropy. 

We have conducted several tests of azimuthal anisotropy with our travel time data. We 
performed least squares inversions to detennine 28 and 40 patterns which provide best fits to the 
age-corrected SS-S residuals. The anisotropy indicated by our regression experiments differs 
significantly from the preferred model of Kuo e, al. (19871 bod. in magnitude and in phase (Figure 
5). Our results indicate that for the 28 model the slow direction for SS-S is N4'W and the peak-to- 
peak magnitude of the effect is less than 1 s; for the 48 model the slow directions are N32‘W and 
N58-E and the magnitude is 2.5 s; for the join, 20 and 48 model tire slow direction is N32'W and 
the magnitude is jus, under 3 s. Km el al. [1987] obtained a peak-to-peak variation with azimuth of 
5-7 s and a slow direction at N13°W. The slowest residuals in the Km el al. [1987] study were 


from north-south paths, i.e., 


nearly along the ridge, and the fastest residuals were from northeast- 


southwest-trending padre with bounce points norih of tire Azores-Gibraltar plate boundary (an area 
noted to be anomalously fas, in their study), so their reported anisonopy may have been at least 
panly the result of unmodelled upper mantle heterogeneity. Our invention for a 28 pattern of 
anisotropy provided a variance reduction of only 2%, compared with 20% for a 48 pattern, and 
22% for a combined 20 and 48 pattern. On the basis of these values of variance reduction and tire 
number of free parameters involved, our results suggest that there is no single coherent pattern of 
upper mande anisonopy in the norih Atlantic. The latest anisotropic upper mantle models obtarned 
horn surface wave tomography [Monlagner and Tanimolo, 1990, also show a complex pattern of 
amsotropy in die region. Any azimuthal anisonopy in tire asthenosphere induced by plate motions 
in tire norih Atlantic may be heterogeneous trecause tire three plates in tire region are slow-moving 
and the return flow is not closely related re plate divergence ( Hager and O'Connell, 1979, 1981; 

Parmentier and Oliver, 1979]. 
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Spatial Patterns of Age-corrected Residuals 

After removal of the dependence on seafloor age, a plot of SS-S navel dme residuals a. the 
SS surface reflection point (Figure 6) shows several interns, ing features. Perhaps the most smiting 
is that residuals in the western Adantic norfh of about 35’ N am on average nearly 4 s mom neganve 
than those to the south. This feature is also nodceable in Figure 3 but is more obvious after age- 
dependence is removed. A similar change a. approximamly this latitude was noted for SS-S 
residuals with bounce positions on the eastern side of tire Mid-Adandc Ridge by Kuo oaf. [19871 
and was attributed to a change in upper rmntie smreture across the Azores-Gibraltar plate boundary. 
The signal we observe is predominantly from data with bounce points on die western side of the 
ridge A map view of the azimuthal distribution is shown in Figure 7 and serves as an aid to assess 
qualitively the geomeuy of wave paths to the south and north of 3TN. We examined the possibiHty 
tha, this signal may be from the Caribbean anomaly, a region of anomalously high velocity in tire 
mantle between 600 and 1400 km depth beneath the Caribbean originally reponed by Jordan and 
Lynn [1974] and further confirmed by Grand [ 1987). If the first leg of tire SS rays propagating ro 
western Europe were to bonom in tire high velocity Caribbean region, the result would be early SS- 
S residuals. This would produce a feature of opposite sign from that observed, so we drscount rt as 
an influence here. Another possible explanation for tire long-wavelength signal could be azimuthal 
anisotropy, but the examination above of possible patterns of azimuthal anisotropy does no, support 


this suggestion. 

Another distinctive feature of tire residuals in Figure 6 is a row of negative values which 
trends nortirwes, to southeast along tire trend of tire New England Seamounts and across tire ridge 
,o the vicinity of the Great Meteor Seamount. This feature comes from even, -station pairs al a 
number of different azimuths and distances so cannot be attributed to a source or receiver effect. 
We do no, observe distinctive anomalies in tire vicinity of tire Bermuda, Azores, Iceland or Canary 
Islands hotspots. The data density is poor for the Bermuda region and Iceland, however, and any 
signal associated with tire Canary Islands may be obscured by the ocean-continent transition. 
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Recently active hotspot islands might be expected to display strong positive (late) residuals, such as 
Stewart and Keen [1978] observed for PP-P residuals at the Fogo Seamounts. In contrast, 
Woodward and Masters [1991] found mostly negative (early) SS-S residuals in the vicinity of the 
Hawaiin hotspot, and Jordan [1979] and Sipkin and Jordan [ 1980] have suggested that the net 
effect of hotspots may be to produce early arrivals because of the presence of high velocities in a 

depleted mantle residuum. 

There is a systematic variation of SS-S residual with latitude, i.e., effectively along the 
direction of the Mid-Atlantic Ridge axis. Age-corrected SS-S residuals with SS bounce points on 
lithosphere younger than 100 My are shown versus latitude in Figure 8. The along-axis variations 
show a variety of scales, notably at wavelengths of about 1000 - 2000 km in the region from 15* to 
35*N, and at about 6000 km wavelength from late (positive residuals) in the south (20-35*) to early 
(negative residuals) farther north (45-55*N). The largest of these variations are robust with respect 
to selective removal of portions of the data. The Iceland region appears as a local maximum 
(positive SS-S delay) on the profile, but the Azores hot spot does not have a distinct seismic signal. 

JOINT INVERSION OF TRAVEL TIME RESIDUALS AND GEOID AND DEPTH ANOMALIES 

Long-wavelength variations in shear wave velocity of the sort depicted in Figure 8 
presumably are a consequence of some combination of variations in temperature and composition of 
the upper mantle. Such lateral variations should also have signatures in other physical quantities 
measurable at these wavelengths, notably gravity (or geoid height) and topography (or residual 
bathymetry), because of the dependence of these quantities on bulk density. Travel time residuals, 
geoid anomalies, and residual depth anomalies are independent quantities dependent in different 
ways on temperature, bulk composition, and their variation with depth. We therefore seek a 
quantitative procedure for treating travel time residuals jointly with geoid and bathymetry data and in 
particular for a combined inversion of all three quantities for horizontal variations in upper mantle 
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temperature and composition. 

To ensure complentarity of data sets, bathymetry and geoid height values are obtained at each 
SS bounce point, and both are corrected for subsidence with seafloor age by means c. the plate 
cooling model [Parsons and Sclater, 1977; Parsons and Richter, 1980]. In this manner we 
effectively normalize all observations to zero age. Bathymetric data arc obtained from the corrected 
Digital Bathymetric Data Base (5’ grid) [US. Naval Oceanographic Office, 1985]. Geoid data arc 
taken from a combined set of Seasat and GEOS3 altimeter data [Marsh el al., 1986]. Data north of 
70° N were not included in the Marsh et al. [1986] data set due to the high probability of being over 
sea ice, so our analysis below is confined to latitudes less than 70*N. We find that the correlation 
of SS-S residuals with the low order geoid is negative, but that at high order the correlation is 
positive (Figure 9). This relationship may indicate a depth dependence of contributions to geoid 
and travel time (e.g., the long wavelength signal may be a lower mantle effect). Low degree 
harmonics are likely linked to deep-seated density heterogeneities and subducting slabs [Hager, 
1984; Hager et al., 1985]. Since we are interested in upper mantle processes, we filter out the long- 
wavelength component of the geoid by subtracting a reference field [Lerch et al., 1979] expanded in 
spherical harmonics to degree and order 7 and tapered to degree and order 1 1. To provide a 
comparable bathymetric data set, bathymetry is high-pass filtered (comer at 4000 km, cutoff at 6000 
km) to remove long-wavelength trends. Along-axis profiles are constructed from the age-corrected 
and filtered geoid and bathymetry data. 

Profiles of age-corrected travel time residuals, geoid, and bathymetry are compared in Figure 
10. While qualitative correlations among profiles are apparent, we seek to quantify possible models 
of temperature and compositional variations that can match these observations. Oceanic bathymetry 
and geoid height are both sensitive to variations in mantle density at depth. Such variations can be 
either thermal or compositional in origin and, like seismic velocity, are presumably related to mantle 
convection and differentiation. For a given density change, the seismic signature of thermal and 
compositional heterogeneity are of opposite sign, so travel time residuals constitute key information 
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for distinguishing between mechanisms of heterogeneity. 

Inversion for Thermal Structure 

We seek to formulate an inversion for the distribution of temperature anomalies T(x,z) 

(where x is along-axis and z is depth) that can produce the along-axis geoid, bathymetry, and travel- 
time anomalies shown in Figure 10. Topography and geoid kernels were calculated for prescribed 
models of viscosity for an incompressible, self-gravitating, Newtonian mande with free slip at the 
surface and the core-mande boundary. The convecting region is assumed to be overlain by a high- 
viscosity layer 40 km thick. We performed calculations both for a mantle of constant viscosity and 
for a mande with a shallow low-viscosity layer. Topography and geoid anomalies depend on the 
viscosity structure, but the predicted travel times do not. Kernels were calculated using a method 
similar to that of Richards and Hager [1984] except that the solution was direcdy integrated across 
the layers instead of being obtained via propagator matrices [McNutt and Judge, 1990]. 

The inversion is best conducted in the horizontal wavenumber domain. The thermal 
anomalies AT(k,z) at depth are related to the predicted dynamic topography h(k) for wavenumber k 
via an integral of the form 

f Zmax 

Ah (k) = J^L- H (k,z) AT (k,z) dz ( 3 ) 

eo-p* L 

[Parsons and Daly, 1983] where a is the volumetric coefficient of thermal expansion, po and p w 
are the densities of the mande and of water at standard temperature and pressure, and z min and z ma x 
are the upper and lower boundaries of the layer in which temperatures are allowed to vary. Table 2 
contains a summary of the constants adopted here. The depth and wavenumber-dependent 
topography kernel H(k,z) is calculated from the equations of continuity and motion given a set of 
boundary conditions, a viscosity model, and a constitutive relation between stress and strain 
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[Parsons and Daly, 1983]. Similarly, the kernel G(k,z) for the geoid relates the thermal anomalies 
to the geoid N(k) via 


AN(k) = 27lF ^ — J G(k,z) AT(k,z) dz 


Z mi x 


gk 


(4) 


[Parsons and Daly, 1983] where T is the gravitational constant, and g is the surface gravitational 
acceleration. 

Sample geoid and topography kernels calculated for different wavenumbers and viscosity 
structures are shown in Figures 1 1 and 12. Cartesian kernels are used throughout this study 
because of their computational efficiency and straightforward application to Fourier transform 
techniques. We have compared extrema of the upper mantle portions of the geoid and topography 
kernels for a layered cartesian Earth and a spherical Earth for a number of wavelengths and different 
viscosity structures (Figure 12), and we note good agreement even at very long wavelengths 
(spherical harmonic order (= 6). This agreement suggests that the results presented here should be 
applicable to the spherical Earth without introducing unreasonably large errors. 

Temperature perturbations at depth can be converted to a seismic velocity perturbation by 
assuming a value for the partial derivative of shear wave velocity with respect to temperature, 
8vs/9T. The resulting two-wave travel time perturbation is given by 


At 


„ dv s f 

I 


^max 


AT (k, z) dz 


v s (z) (1 - p 2 v s (z) ) 


2.1/2 


(5) 


where v s (z) is from the reference shear velocity model [Dziewonski and Anderson, 1981] and p is 
the ray parameter, generally taken to be a representative value for the range of epicentral distances 
considered here. We use a value of -0.6 m/s K'l for dv^T. This value is higher than the values 
of Anderson et al. [1968] and Kumazawa and Anderson [1969] at standard temperature and 
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pressure but is similar to the value of -0.62 m/s K’> determined by McNutt and Judge (19901 by a 
least squares fit of Love-wave phase velocities to predicted temperature of the lithosphere. Such a 
value is consistent with the change in P-wave velocity with temperature, dvpOT = -0.5 m/s K >, 
found from modeling wave propagation along subducting slabs [Creager and Jordan, 1986, Fische 
e,al„ 1988] if we assume that 3vs/3T = 1.1 dvpOT [Woodhouse and DziewonsU, 1984]. Partial 
melt would increase the value of tivsOT [Sleep, 1974; Sato and Sacks. 1989], but simultaneous 
analysis of both shear and congressional differential travel times by Woodward and Masters 1 1991 ] 
indicates that significant pamal melting is no. required to explain the differential travel time residuals 

in the north Atlantic region. 

The forward problem consists of calculating geoid, topography, and travel ume residual 
profiles given a starting twodimensional temperature structure T(x,z). The inverse problem 
consists of finding a temperature structure that predicts (via equations 3-5) geoid, topography, and 
navel time profiles which best fit those observed. The familiar matrix equation d = A m is formed 
from discrete versions of equations 3 - 5. The data vector d consists of the topography, geoid, and 
navel time residuals, the model vector m contains die temperature variations for which we are 
solving, and the matrix A contains the coefficients and kernels which relate the data to the model. 

As a check on our procedure, we constructed a forward problem for geoid and topography and 
found good agreement with the modelling results of McKenzie etal. (1980], 

The bathymetry, geoid, and navel time profiles of Figure 10 are interpolated to a constant 
spacing, demeaned, tapered a. bod. ends with a 10% sine squared taper, and Fourier transformed. 
Since our profile extends from 10 to 72»N, tire firs, and las. 10% of the profile (10 - 16°N and 66 - 
72°N) will be affected by the taper. The 3n x 1 data vector d is then constructed, using tire complex 
(to rerain both amplitude and phase) bathymetry, geoid, and navel time data sampled a. n discrete 

wavenumbers: 

d = [ Ah(ki),...,Ah(k n ),AN(ki),...,AN(k n ),At(ki),..., At(k n ) ] T ^ 
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where T denotes transpose. For the case 

single layer, the n x 1 model vector m is given by 


where temperature perturbations are constrained to be in a 


m 


= [ AT(ki),...,AT(k n )l T 


( 7 ) 


For the more general case of a 


multi-layer system, the nj x 1 model vector m is given by 


( 8 ) 


m . t Am, x,) ATfki.Zj). AT(kazt) AT(k 2 z j) .AT<k,z„...,AT(k n Zj)]? 

where z, is the layer index and j is the .oral number of layers. In this paper we perform tnverstons 
for single-layer models only. The “layers” of temperature variations are tndependen, of the 
“layering” system of lid, low viscosity zone, and mantle which we use for the calculation of 
kernels, although major changes in viscosity would rend to segment AT as well. The temperature 
layering simply refers to that region bounded by zmin and z max in die integrals of equations - . 

The 3n x n matrix A contains the coefftcients and kernels that relate the temperature 
perturbations to the observations, which for the single-layer case is given by 
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P 00 - £ H(ki,z) ta 0 

po ■ P W wmam 

0 p< ^ - £ HOtn) Az 

pO - P* Mail 




(9) 


m m aui* A contains both bathymetry and topography kernels and is Utus viscosity dependent; 
u „ a viscosity smteture must be assnmed. We solve the equation d = A n, by leas, squares 


m = ( A R dd A ) A R dd d 


( 10 ) 


where A is the complex conjugate transpose of the matrix 
matrix Rdd is discussed in Appendix B. Equation (10) is 
variance reduction is calculated via 


A. Construction of the data covariance 
solved for the solution vector m, and 
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variance reduction = 1 - 


(d-Am) Rdd (d-Am) 
d Rid d 


(ID 


The resulting model vector m is inverse Fourier transformed back to the spatial domain to produce 
an along-axis temperature profile. The solution m is also substituted into equations 3 - 5 to 
compare predicted geoid, bathymetry, and travel time residuals with those observed. 

Six inversion experiments for temperature structure were performed (Table 3). Inversions 
were carried out for two different viscosity models and for three different thicknesses of the layer in 
which lateral temperature variations were assumed to occur. Because topography and geoid 
anomalies depend only on the ratios of viscosity in different layers [Robinson et al . , 1987; Hong et 
al., 1990], we set the dimensionless viscosity of the layer representing the bulk of the mantle to 
unity. In one viscosity model, termed the “constant viscosity mantle,” a 40-km-thick high-viscosity 
lid overlies a unit viscosity mantle. We set the viscosity in the lid to 10 4 , which effectively mimics 
rigid behavior. In a second model, a 160-km-thick low-viscosity zone is present beneath a 40-km- 
thick lid; the viscosity in the low-viscosity zone is a factor of 100 less than in the underlying mantle. 
The thickness of the layer of temperature perturbations was taken variously to extend from 0- 1 50 
km depth, 0-300 km depth, and 0-650 km depth. The matrix A is different for each of these cases, 
as it involves viscosity-dependent geoid and topography kernels and also a summation over depth. 

Inversion results for the constant-viscosity-mantle cases are shown m Figure 13. The 
“observed” profile is actually a filtered version of the observations, containing only the wavelengths 
used in the inversion (1400 to 7 100 km). Predicted profiles were calculated from equation 5. For 
these solutions, the long-wavelength fit to geoid is better than at short wavelengths. The fit to 
bathymetry is poor. The predicted magnitude of the SS-S residuals range from a factor of 5 too 
small for the 650-km-thick layer to a factor of about 1.5 too small for the 150-km-thick layer. 
Increasing the temperature variations to improve the fit to the SS-S residuals leads to predicted 
geoid variations that are too large. The highest total variance reduction and best fit for the constant- 
viscosity cases come when lateral variations are constrained to shallow (0-150 km) depth. The 
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variance reduction is 25% for bathymetry, 79% for geoid, and 58% for travel times. The total 
variance reduction is 53%. The variation in temperature is 180 K for the 150-km-thick layer, and 

only 60 K for the 300-km-thick layer. 

Figure 14 shows inversion results for the models with a thin low-viscosity zone. A good fit 
lo both geoid and travel time is found, although the alignment in phase of predicted and observed 
geoid is not as good as for the constant-viscosity case. The fit to bathymetry is again poor. The 
total variance reduction for the 150-km-thick and 300-km-lhick layers arc both 57%, although the 
shallow model provides slightly higher variance reduction for bathymetry (27% for 0-150 km deep 
layer, 24% for 0-300 km deep layer) and the 300-km-thick-layer model provides higher variance 
reduction for geoid (79% for 0- 150 km deep layer, 85% for 0-300 km deep layer). The variation in 
temperature for the 150-km-thick layer is 230 K and in the 300-km-thick layer is 1 10 K. 

We have explored the hypothesis that the lack of correlation of predicted and observed 
topography is an indication that the source of variations in the geoid and travel time anomalies is 
deep. To test this hypothesis, we performed inversions with temperature variations restricted to 
deeper layers and found that fits to topography were still poor. Il is possible that the bathymetric 
signal is dominated by crustal thickness variations which are not included in our calculation of dy- 
namic topography. An assessment of such thickness variations is discussed further in Appendix B 


Inversion for Compositional Variations 

A possible alternative to along-axis variation in mantle temperature is lateral variation in bulk 
mantle composition, due perhaps to a variable extent of melt extraction or different degrees of 
mixing of compositionally distinct volumes of mantle material. The dynamical effects of 
compositional* induced density variations can be large [O’Hara, 1975; Boyd and McCallister, 
1976; Oxburgh and Parmentier, 1977; Sotin and Parmentier, 1989]. The fraction of mantle 
potentially extractable as basaltic melt is thought to be 15-25% [e.g.. Green and Uebermann, 
1976]. Thus, for every volume of basalt removed from the mantle, a volume of residuum several 
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dmes larger is left behind. The effect of basalt depledon is to increase the molar ratio Mg/(Mg + Fe) 
(or Mg#) in the residuum, which reduces the density and increases the seismic velocities [e.g., 
Uebermann, 1970; Akimoto, 1972). For example, subtraction of 20 mole % olivine basalt from 
pyrolite can decrease the density of the residuum by neariy 2%. equivalent to a thermal perturbation 
of nearly 500 K [ Jordan , 1979). Thus compositional changes need only be slight to produce effects 
on the older of 100 K, comparable to values obtained from the inversions for temperature 
variations. In this section we explore the effects of compositional variations parameterized in terms 
of the variation in the Mg# in the upper mande along the ridge. Our motivation for parameterizing 
compositional variations simply in terms of Mg# is that differences in this quantity yield significant 
variations in seismic velocity and density, in conirast to most other measures of degree of melt 


extraction. 

Partial derivatives of density and seismic velocity with respect to Mg# are obtained from 
Akimoto [1972]. These values were measured on a suite of samples ranging from pure forsterite 
(MgzSiOa) to pure fayalite (Fe 2 Si0 4 ). While these partial derivatives are at standard temperature 
and pressure, it is expected that a change to elevated temperature and pressure wUl have only a 
second order effect, since temperature and pressure conections work in opposite directions [Jordan, 
1979). Above the solidus temperature, however, the amount and distribution of partial melt, which 
may depend strongly on composition and particularly volatile content, is important The presence 
of melt is likely to have a larger effect on shear wave velocities than on bulk density. Calculations 
of melt migration, however, suggest that once created, melt segregates rapidly by a percolation 
mechanism [e.g. Scott and Stevenson, 1989], so that the melt fraction present in the mantle at any 
given time is probably small. Studies of mantle peridotites [Johnson et al, 1990] also support the 

importance of fractional melting. 

I. is straightforward to conven equations (3) and (4) to relations between geoid or 
topography and a compositionally induced density perturbation by means of the relation 
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Ap = - p 0 a AT (12) 

Compositional anomalies at depth yield a dynamic topography h(k) given by 


^max 

Ah(k) = - 1 — f H(k,z) AMg(k,z) dz (13) 

P 0 “Pw oMg J 

Z min 

where AMg represents the fractional change in the Mg#. Compositional anomalies yield a geoid 
anomaly 


Z max 

AN(k) = ^f4dg j ° (k,Z) AMg ^ k,Z ^ dZ 


(14) 


For a compositional perturbation at depth the resulting two-wave travel time perturbation is given by 


Zimx 


At(k) = 2 


3v s 


J 


AMg(k,z) dz 


3Mg J v ( z ) 2 (i . p 2 v s (z) ) 


2.1/2 


(15) 


Using equations (13) - (15), an inversion scheme similar to that used for thermal 
perturbations is formed. The solution vector now has the form 

m = [AMg(ki) AMg(k n )] T (16) 

The data vector remains the same as in equation (6), while the matrix of coefficients, A, changes to 
reflect the relation between the data and mantle composition, rather than temperature, as outlined in 

equations (13) - (15). 

The results of the inversions for compositional variations are summarized in Table 3 and in 
Figures 15 and 16. We are unable to match simultaneously both SS-S travel time residuals and 
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geoid and bathymetric anomalies with so.e.y mantle compositional variations for either a constant 
viscosity mantle or one with a low viscosity zone. This is no, apprising, as the travel times are for 
the most part positively related with geoid and bathymetry. bu, composidonal vanattons (a, leas, 
for the MgaSiOa - FeaSiOa system examined hem) have an opposite effort on travel dme and geotd- 

bathymetry. 

for ^ constant viscosity case, the fit to the geoid is excellent, and the fi, to bathymetry ts 
slightly better than in die invention for temperature. The fit to SS-S residuals is so poor that die 
variance reduce is negative for travel dme. Urge compositional changes would be tequued to 
affect travel times, whereas only small compositional changes am needed to ptuduce signtfican, 
density contrasts to match the geoid signal. The total variance reduction for die constant viscostty 
case does no, vary gready (ftom 32-33%) for compositional changes constrained to be over 
diffemn, depth intervals, though the variance reductions for individual da. sets (badiymeuy, geoid, 
travel time) vary significandy from mode, to model (see Table 3). The range in Mg# is about l%if 
die variation is constrained to die depth range 0-150 km and only 0.1% for die 0-650 km depth 

range. 

Figure 16 shows inversion results for the model with a low viscosity zone. A good fit to 
tod, geoid and bathymetry is found, although the alignment of predicted and observed geotd ,s not 
- - 600 d as in die constant viscosity case. The fi, to bathymetry is die best of any models so far. 
The total variance reduction is still low (43 to 49%), due to die poor fi. til navel times (negattve 
variance reduction in all cases except die 0-650 km model). The range in Mg# is 2.4% if 
constrained to 0-150 km depth, 1.3% over 0-300 km depth, and 0.5% over 0-650 km depth. 


Joint Inversion for Temperature and Composition 

We next explore whether a combination of tempereture and compositional variations can 

pmv.de a good match to die observed geoid, nave, time, and b,hymeny. loin, inversions provide 

improved fits to all da. a, die expense of inducing additional free parameters. For these 
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inversions the to vector remains the same as in equation (6). the solution vector is modified to 
include botit temperature and composition, and the matrix of coefficients. A, includes the effects of 
both temperature and compos, tion. The matrix-building equations become, for example, for 

topography. 


Zmax 

f“ x AT(V d7 + \ — L f H(k,z) AMg(k,z) dz (17) 

_ H(k,z) AT(k,z) dz + p Pw 3Mg J 

Y J u z . 

Z min 

min 


Po a 

which is simply a combination of equations (3) and (13), The new geoid equation comes from a 
combination of equations (4) and (14) and the travel time equation from a combination of equations 
(5) and (15). Cross terms, such as compositional changes induced by increases or decreases 

temperature, are neglected. 

The results for the join, inversion for temperature and composition are summarized in Table 3 
and Figures 16 and 17. The travel time residuals are well-modeled in all cases, as are the geoid 
to. The topography is best fit for the case with a low viscosity zone. Resolution of the depth 
interval of tire most impomnt lateral variations is rather poor. The topography is fit manually 
betier for the case where temperature and compositional anomalies are constrained to be shallower 
than 300 km. For the constant viscosity mantle, the temperature variations range from 210 K, if 
constrained to 0-150 km depth, ,o55K if over 0-650 km depth; variations in Mg# range from 1.5% 
if over 0- 1 50 km depth to 0.4% for 0-650 km depth. For tire case with an upper mantle low 
viscosity zone, tire temperature variations are similar to those in tire constant viscosity case, bu, the 
variations in Mg# are larger, from over 2% for 0-150 km depth to nearly 1% for 0-650 km depth. 

The navel time residuals are perfectly fit in the join, inventions for temperature and 
composition (Table 3). This occurs because of the way tire model parametera act in a similar 
manner on both geoid and bathymeny, producing a singular matrix if only geoid and bathymetry 
to are inverted for both temperature and composition. Undamped leas, squares always pmvides 
perfect solutions when the number of equations is equal to the number of unknowns unless the 
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matrix to be inverted is singular. If we perform an inversion including only travel time and geoid 
data, we have the same number of equations as unknowns, the matrix is nonsingular, and we obtain 
perfect fits to both travel time and geoid Similarly, if we perform an inversion of travel time and 
bathymetry data, we again obtain perfect fits to both data sets. If we perform an inversion of geoid 
and bathymetry data, however, we are unable to obtain solutions without applying damping. In the 
joint inversion of travel time, geoid, and bathymetry, we have more equations than unknowns and 
the inversion is overdetermined. However, the travel times are perfeedy determined in this case 
because of the nonuniqueness inhetent with geoid and bathymetry. We have performed undamped 
inversions with various weightings on the geoid, bathymeny. and travel time data, and in all cases 


the travel times remain perfectly fit. 

We have also performed joint inversion for temperature and composition with Mg# variations 
constrained to be in the upper 50 km of the lithosphere so as to mimic compositional variations due 
solely to variable melt extraction at the ridge. Temperature perturbations were allowed to remain 
wilhin the depth ranges adopted earlier. The results for this invention are summarized in Table 3 
and Figures 19 and 20. Tbe variance reduction was similar for the constant viscosity case and for 
the model with a low viscosi ty zone. In general, the geoid is fit very well, the predicted amplitudes 
are a bit low for travel time residuals, and the topography fit is slightly out of phase. For the 
constant viscosity mantle, the range in temperature is 210 K over 0-150 km depth and 25 K over 0- 
650 km, while Mg# variations constrained to be confined to 0-50 km depth were over 5%. For the 
case with a low viscosity zone, the temperature variations were not dramatically different from those 
in the constant viscosity case, and variations in Mg# were about 4.5%. The inversion solution 
shows high temperatures near 30'N and low temperatures in the region from 50-60’N. Iceland also 
appears to be underlain by high-temperature mantle. Going from south to north along tile ridge, 
compositional variations indicate low Mg# in die vicinity of 20-30‘N. high Mg# in die Azores 
region (40*N), low values near 50*N, and high values near 60 N. 
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DISCUSSION 

The temperature and compositional variations in Figures 13-20 are broadly consistent with 
observed travel time, geoid, and bathymetry anomalies in the norih Atlantic region. Temperature 
variations alone can account for most of the observed anomalies. In contrast, compositional 
variations alone cannot match all anomalies simultaneously. We infer that a component of the 
observed anomalies is due to long-wavelength variations in upper mantle temperature. Joint 
inversions for temperature and composition provide better fits than single-variable models, but at 

the expense of introducing additional free parameters. 

It is difficult to select a ‘best' model from the suite of inventions presented. The variance 

reductions in Table 3 serve as a guide, but independent criteria may allow us to reject some of the 
models, even those with high variance reductions. In particular, those models with large 
temperature variations (well in excess of 100 K) can be seriously questioned. Lateral temperature 
variations at upper mantle levels beneath oceanic ridges are thought to be no more than about 300 K 
globally [ Klein and Langmuir, 1987; Johnson et al„ 1990], so a variation in temperature of 230 K 
(as in the inversion with a low-viscosity aone and a 150-km-thick layer of temperature 
perturbations) solely within a section of the north Atlantic is probably unreasonably large. Further, 
as While and McKenzie [ 1989] have noted, relatively small increases in mantle temperature above 
values typical for the mid-ocean ridge are sufficient to cause large increases in melt production. 
Their models indicate that, for fixed bulk composition, an increase of 100 K above normal doubles 
die amount of melt while a 200 K increase can quadruple it. Such increased melt production should 
lead to approximately corresponding increases in ctustal thickness. Variations in oceanic crustal 
thickness away from fracture rones, however, are generally thought to be small, with thicknesses 
typically 6-7 km and ranging from 4.5 to 8.5 km [ Spudich and Orcutt, 1980; White, 1984; Purdy 
and Detrick, 1986]. In the joint inversion for temperature and composition, temperature vananons 
if confined to 150 km depth are excessive (over 200 K) and if the variations extend over 0-650 km 
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the fit to topography is poor, especially for the constant-viscosity mantle. On the basis of these 
results we prefer the models with temperature variations occutring over 0-300 km depth. For the 
constant viscosity mantle, the temperature variadon is 1 10 K, and the variation in Mg# is 0,75%. 

For the case with an upper rode low viscosity zone, the predicted temperature variation is 125 K, 
and the variation in Mg# is 1.1 ft. the ratal variance reduction is greater in the model with a low 

viscosity zone. 

Even a temperature variation of about 100 K is high for a mantle of constant composition, 
since we do not observe increased crustal thickness in regions that our models indicate have high 
temperatures. The assumption of app.oxima.ely constant upper mantle composition warrants 
discussion. In particular, lateral vanation in trace amounts of mantle volatiles may have a large 
effect on seismic velocity at a given temperature. The presence of even a slight amount of water, 
for instance, is sufficient to cause a significant decrease in the initial melting temperature of 
peridotite [Wyllie, 1971). Estimates of volatile contents and their lateral variations in the north 
Atlantic region have been made from measurements of abundances of halogens, SiOz, K 2 0, and 
H 2 o in basalts and from the volumes of vesicles in basalts [Schilling el al„ 1980, 1983; Schilling. 
1986; Michael, 1988]. These studies indicate that Cl, Br, F, and H 2 G contents increase toward the 
Azores and Iceland and that H 2 0 is two to three times more abundant in Mid-Atlantic Ridge basalts 
erupted over the Azores platfotm than a. adjacent notmal ridge segments. The effect of volatiles on 
density and shear wave velocity will be slight a. subsolidus temperatures but can be major over the 
melting interval [Goelze, 1977], The presence of melt will act to decrease significantly the seismic 
velocity [Duschenes and Solomon. 1973] and, to a lesser extent, lower the density of the mantle. 
To the extent that seismic velocity depends on proximity of the temperature to the solidus 
temperature [Sum e, al.. 1988, 1989], volatile content can trade off with temperature in its effect on 
velocity at subsolidus conditions. Thus, variation in volatile content could lessen the variations in 

melt production implied by the inversion solutions. 

Even without significant variations in volatile content, it is clearly an oversimplification to 
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parameterize mantle composition in ternis of only a single quantity. Further we have assumed that 
the partial derivatives of bulk density and seismic velocity with respect to Mg# that are those for 
olivine [Akhnom, 1972], The work of Jordan [1979] indicates drat these derivatives remain nearly 
constant for many different mande compositions (i.e., pyrolite-type composmons with various 
amounts of olivine, orthopyroxene, clinopyroxene, spinel, and garnet), so the latter assumpnon is 
sound. However, at any given Mg#, orthopyroxene and clinopyroxene have lower velocities and 
are less dense than olivine, while game, and spine, are seismically faster and denser than olivine 
[Jordan, 1979], so an increase in the weigh, pereent of orihopyroxene and clinopyroxene or a 
decrease in die weigh, percent of game, and spinel with respect to olivine in die mantle could 
counteract some of die temperature variations obtained under die assumption of effectively uniform 
mineralogy. Several studies [Wood, 1979; Jacques and Green, 1980; Dick et al„ 1984] have 
suggested that compositional variations in die mande are plausible. Indeed a number of workers 
[e.g„ Davies, 1984; AUegre e, al„ 1984] favor dynamic models for the mande in which dispersed 
heterogeneities of various sizes and shapes are passively embedded in a continually mixed, 
convecting mande. Variations in modal fractions of olivine, orihopyroxene, and clinopyroxene in 
peridotites recovered along die Mid-Adantic Ridge have been reponed in several studies [ Dick e.al„ 
1984; Michael and Bonash, 1985], These variations are typically attributed to different degrees of 
mel, extraction but could also be partially due to intrinsic upper mande heterogeneity. For example, 
the relative fractions of olivine, clinopyroxene, and orihopyroxene indicated by Michael and Bonam 
[19851 at 26*N and 30*N, if extended to depth, could counteract a portion of the temperature 

differences indicated by the inversion solutions for these regions. 

Chemical analysis of dredged peridotites in the norih Adamic indicate a range of about 2.5% 
variation in Mg# [Michael and Bonam, 1985], This value is intermediate between wha, we find for 
models with compositional variations constrained to be shallow (4.5 to 6% variation) and those 
models with compositional variations in the same depth ranges as die thermal variations (1-2%), 
This suggests that compositional variations may be concentrated slighdy shallower than die 
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temperature variations. Michael and Bonaai [1985] present an along-axis prof.le of Mg# variations 
tan dredged peridoti.es which can be compared with our emulated profile. The main feature in 
their profile is a zone of high values of Mg# in the Azores region, from 34-45'N, relative to the res. 
of the ridge, consistent with our modelling results. Their data sampling is too sparse to delineate 
other long-wavelength features. Their average value for 26’N also has a high Mg# relative to 
adjacent data. This is consistent with our observation of early SS-S navel times and low geotd tn 
this region. This anomaly is of too short a wavelength « 1000 km), however, to resolve in our 
inversions. We should note that comparisons ment caution, as small scale features, such as t 
due to ridge segmentation, can produce large differences in composition between peridotites over 
scales of tens of kilometers. In addition, dredged peridotites are mostly from fracture zone 
environments, which may not be representative of typical ridge mande [Dick, 1989). 

On the SS-S residual profile the Iceland region appears as a local maximum (late SS) but the 
Azores hotspot does not show a distinct seismic signal. The inversion results for these two regions 
are also markedly different The results of the join, .aversion for temperature and composition 
predict a high Mg# in the Azores region while indicated temperatures are no. anomalously high. At 
Iceland, in contiast, high temperatures dominate. Work by Schilling [1986] and Bonaai [1990] 
oudines the differences in geochemical signatures between the Azores and Iceland ho. spots. These 
workers suggest that Iceland is a "traditional" plume hot spot, with a predominandy .hemal origin, 
but that the Azores might be more aptly named a “we. spot" because of the presence of excess 
hydrous phases and the lack ofathemal anomaly. Bonaai [1990] suggests that because the 
Azores hotspot is rich in volatiles, enhanced melt production could occur with li.de or no increase tn 
temperature. The high Mg# indicated in our inversions allows the region to be seismically fas, (as 
we observe) bu, of low density (as geoid and bathymetry require). The results are consistent with 
the hypothesis that die Azores hot spot is no. associated with a plume-tike thermal anomaly. 
Inversion of surface wave dispersion data can potentially provide furdier tests of these ideas, but 
studies to date have yielded apparendy conflicting results. Results of several such investigations 
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[Nakaniski and Anderson, 1984 -.Tammon, 1990; Ztang tmd ranimorr,, 1990] suggestthatthe 
AzDrc s region is seismical.y slow a. depths .ess than 300 km b„, a study utilizing 50 - 20 (Vs-penod 
Rayleigh waves by Mac**,', al. (.989, dees not. These differs may be pardaUy atmbutable 
,o ure differences in wave periods employed and nrede of analysis fmm study to study. 1. may be 
possible that what appear to be low velocities at die Azores are a result of horizontal smoothing of 
*. low velocities along the ridge and have Utile to do with the actital structine in the Azores region. 
None of these iong-period surface wave studies resolve a distinctive anomaly a. Iceland. Clearly, 
more work is needed to resolve the upper mantle velocity structure of hot spot regions. 

Bonatti [1990] has constructed pmfiles of the equilibrium temperature of dredged pen otites 

•c frr™ n to 60-N bv means of two different geothermometers 
along the Mid- Atlantic Ridge axis from 0 to 00 in Dy meai 

[Wells 1977; Undsley, 1983,. Comparison of these pmft.es with the along-axis temperature 
variations obtained from our inversions reveals a number of qualitative conelations as well as a few 
discrepancies. The range of temperature variations in the ptoftle based on the M U** 
geothermometer is about 150 K, neglecting high values termed -anomalous" When the htgh values 
are included the range increases to 400 K. The profile utilizing die Wells [ 1977] geothermometer 
Has a range of 100 K neglecting the anomalous values and 350 K including them. The hrghes. 
temperatures in our inversions are near 30'N (Figures 17-18), a region showing a sl.gh. peak tn 
Bonatti's temperature profile estimated according to Lindsley ]1983, and a very weak rise in the 
profile utilizing die Wells (1977, geothermometer. mere is a small dip in temperature a. (a 
region which we find to be seistidcally fast) in the Un*,ey ,1983, and Wells ,1977, profiles, but 
tire difference may no, be significant considering the enor bars. Ber gm an end Solomon [ 1989, also 
found die upper mantle near 26'N to be seismically fas. fiom an analysis of teleseismic P-wave 
travel time residuals from eafihquakes in this region reconied by a local ocean-botiom seismic 
network. The .owes, temperatures on the profiles of ,1990, are a, 43*. Temperatures 

fmm our inversion solutions are also low in this region, although the Bonaai ,1990, profiles 
indicate an increase in temperature proceeding nonh from 43‘K to 53’N, whereas our results favor 
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conned low temperatures. Pat, of the differ between out .suits and the geochemical sutdies 
may be anributed to the fact that the depth sampledby basalts and peridoti.es » U* » be 
shallower than the layer thicknesses of most of our models. Assuming that the 6-km-thick oceanic 
crust was formed by 9 to 22 % partial melting of the mantle l Klein arui^muir, 1987], then the 
volume of residual peridotite will extend from the base of the crus, tit somewhere between 30 and 
70 km depth. The amount of depletion will vary with depth if we assume a fractional melung 
model. Our models with compos, tional variations confined to depths less titan 50 km are most 

representative of shallow fractionation and differentiation. 

Several improvements in future studies of the type ptesented here may be envisioned. Our 

models thus far have been limited to simply parameterized one-dimensional variations in 
temperature and composition within a single layer. I. is likely that these lateral varrauons are no, 
constant within a given layer and .ha, there are two-dimensional lateral variations .dependent of 
lithospheric aging. The rechniques outlined in this paper can be generalized » a multilayer system 
and to two-dimensional wavenumber (see equation 8). bu, we do no, fee, that the resolution of our 
data can justify more complicated mcxlels a. this time. Kernels for seismic surface waves are 
strongly peaked in the upper mantle, and such da. would provide a useful consent tn future 
models. The inclusion of surface wave da. would help to distinguish between lithosphenc and 
astitenospheric effects and may allow for two or more independently resolved layers. Extension of 
die modelling to three dimensions would permit an assessment of the degree to which mantle 
anomalies beneath the ridge extend off axis. Implicit in our age-conection is the assumption tita, tile 
anomalous pmperties of the ridge mantle are steady s.t= on a time scale of 100 My. Recent 
seafloor surveys and theoretical studies (e.g., Poctolny e, a,.. 1988; Scon Stevenson, 1989, 
bring tins assumption into question and suggest tiia, a, leas, on shore time scales « 1 My) and a, 
slow spreading rates (as in the Atlantic) intermittent periods of melting and crustal forenauon may be 
separated by periods with Utile or no mel, pnMuction. These temporation variations are likely to be 
averaged outi however, over the typical horizontal wavelength (100 km) of a long-period SS wave. 



33 


Another limitation of our models is that they depend on the assumed values of several 
physical constants. It is straightforward, however, to estimate the effect of choosing different 
values. The viscosity structures we employ are also quite simple but have been chosen to represent 
two models widely invoked in other studies - a constant or nearly constant viscosity mantle (e.g., 
Peltier, 1989] and a mantle with a thin low viscosity layer le.g.. Craig and McKenzie, 1986; 
Robinson et al„ 1987]. The viscosity structure of the Earth may be temperature and pressure 
dependent or vtuy laterally, but we have not considered viscosity structures of this type. Nor have 
we modelled the effects of partial melting which could accompany the temperature variations we 
predict The effect of retained melt on the physical properties of the mantle depends critically on the 
melt fraction and geometry, characteristics presently poorly known. Sato et at. [1988, 1989] 
downplay the importance of partial melt and suggest that most mantle seismic velocity anomalies 
can be explained by temperature variations at subsolidus conditions. The combined analysis of both 
shear and compressional differential travel times also suggest that significant partial melting is not 
required to explain the travel time residuals in the north Adamic region [Woodward and Masters, 

1991]. , 


CONCLUSIONS 


We have measured 500 SS-S differential travel times for paths in the north Atlantic region. 
The SS-S travel time residual decreases linearly with square root of age, in general agreement with 
the plate cooling model to an age of 80-100 My [ Parsons and Sclater, 1977]. Azimuthal anisotropy 
is not clearly resolved, and the azimuthal patterns of our data are not consistent with the preferred 
upper mande anisotropy model of Kuo et al. [1987] for the north Adantic. An along-axis profile of 
age-coirected travel time residuals displays significant long-wavelength variations, notably at 
wavelengths of 1000-2000 km. The largest of these variations are robust with respect to selective 

removal of portions of the data. 
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We have formulated a join, inversion of travel time residuals and geoid and bathymetric 
anomalies for lateral variation in upper mande temperature and composition. On the basis of 
variance reduction, inversion for temperature favors the presence of an upper mande'tow viscos.ty 
^ and temperature anomalies concenuated a, depths less than 300 km. We are unable to match 
navel time residuals simultaneously with geoid and bathymetry solely with lateral vanatrons m bulk 
composition (Mg#). Join, inversions for temperature and condition provide good fits to both 
navel time and and geoid regardless of viscosity structure or layer depth and thickness, but the best 
frts to bathymetry come ftom models with a low-viscosi ty none and thermal or compositional 
variations confined to shallow depth. The Mg# variations predicted in the join, inversion for 
temperature and composition are comparable to those found by Michael and Bonam [1985] in a 
study of dredged peridotites along the Mid-Atlantic Ridge and may be related to variations in melt 

production along the ridge. 

IT* preferred inversion solutions have variations in upper mantle temperature along tire Mid- 
Atlantic Ridge of about 100 K. For a constant bulk composition, such a temperature variation 
would produce about a 7 km variation in crustal tirickness [White end McKentie, 1989], iarger man 
is generally observed Vpudich and Orem, 1980, Whi, ,e. 1984; Purdy and Derick, 1986], 
Introducing compositional variations as well as temperature variations in tire inversions does no. 
change the range of temperature appreciably. Tire presence of volatiles in the mande can have a 
strong effect on temperatures required for melting, and variations in volatile content along me ndge 
may reduce me large variation in melt production implied by me lateral temperature variations 

indicated in our models. 
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APPENDIX A: ESTIMATION OF ERRORS FOR SS-S DIFFERENTIAL TRAVEL TIMES 

It is important to quantify the uncertainties in the differential travel time measurements. After 
cross -correlation, the “quality” of each individual SS-S measurement is rated and a grade is 
assigned. The cross correlation coefficient, which describes the degree of fit between the synthetic 
and real SS phases, is used as an objective aid in the assignment of quality. However, our final 
assignment of quality is largely subjective and based upon visual inspection of the "synthetic" SS, 
real SS, and cross correlogram, taking into account the sharpness of the arrivals and their 
alignment, the clarity of the seismogram, and the appearance of a single clear peak in the cross 
correlation function. An "A" quality grade indicates an excellent fit, "B” quality indicates good 
phase alignment but only a fair fit, and a "C" quality grade indicates a poor fit or some ambiguity as 
to phase alignment In addition to A, B, and C grades, there were data that were rejected due to 

poor signal to noise ratio for either the S or SS phases. 

Assuming that the uncertainty in an individual measurement comes from a combination of 
measurement eiror, unmodeled lower mantle structure, and epicentral error, we write, for example, 
for the measurement variance of an "A" quality datum: 


a A 2 = °Am 2 + °lm 2 + G epi 2 

where a A is the total uncertainty, o Am is the measurement error, is the uncertainty due to 
unmodeled lower mantle structure, and o epi is the epicentral error. We assume that and c epi 
are the same for A, B, and C quality measurements, but the measurement error is obviously a 

strong function of data quality. 

Effect of Epicentral Error 

In general, epicentral errors affect ihe travel times only slightly. The events used m this 

C-3 
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study well recorded by a large number of stadons over a wide range of animuths, and typ.cal 
epicentral mislocations am pmbably less dtan ,0 km (which would yield a differential iravel-dme 
errorof 0.35 s a, 75’ distance). The navel dmes are even less sensitive to enors in focal depth; an 
enor in depth of 25 km conntbutes only afctut 0.3 s to the SS-S res, dual. Using tire rule of thumb 
that one standard deviation is about one half of the estimated extremes, we adopt o epl = 0.75 s as a 
conservative estimate of epicentral error. 


Effect of Unmodelled Lower Mantle Heterogeneity 

We estimate the likely magnitude of lateral variations in the shear wave velocity of the lower 

mantle fan models of lower mantle heterogeneity in P wave velocity (such as model L02.56 of 
Dzim onski (1984,). The average variation in navel timesof dime. P waves bonoming in the lower 
mantle is in the range ± 0.5 s. Global tomographic studies by Dziewonsld and Woodhouse [1987] 
indicate that the scaling tatio (SV^ /(SVplV p ) - 2 in the lower mantle. Such a scaling is also 
suggested by comparison of lower mantle P wave mtrfels with the mcen. lower mantle S model of 
Tanimow [1990], Assuming such an S to P velocity anomaly scaling, the resulting variatton m S 
wave arrival time contributed by the lower mantle would likely be about ± 1.5 s, a (taction of the 
observed range in SS-S residual. While the major features of lower mantle model L02.56 
Vronski, 1984] and the lower mantle portions of Tan.mo.oS [1990, model am for the most 
pari similar, enough differences exist drat the application of a lower mantle “correction” to our data 
nugh, add more uncerrainty than i, remove, Furiher, absolute S-wave navel times do no, show 
enough variance for us to suspect large lower mantle effects (e.g., Randall, 1971; Girardin and 
Poapine., 1974; Har, and Bader, 1978; UhrHananer, 1978, 1979], and the work of Gudmundson 
e, al. [1990] indicates that most of the variance from dm ISC rabies is atiriburable to the shallow 
mantle, i.e„ most of the Earth’s heteregeneity is in the upper mande and the lower mantle is fanly 
homogeneous. On the basis of the above information, we set o lm = 0.5 s for our study. 
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Measurement error 

As an objective means to obtain error estimates, we examine the scatter in A, B, and C 
quality picks in a small region. We measured the mot mean square (rms) difference between travel 
time residuals of the same grade (A, B, or C) with bounce points separated by less than 80 km and 
with differences in path azimuth at the bounce point of less titan 10'. An 80-km distance is less 
than the horizontal wavelength of SS (which is about 1 80 km at 25 s period) so we do not expect 
much contribution to the rms difference from actual lateral variations in structure. The mis 
difference for the 16 A quality residual pairs which were within 80 km of each other was 1 . 15 s. 

. uality picks, an rms difference of 2.08 s was measured using 20 residual pairs, and for C 

quality picks 44 residual pairs yielded an rms difference of 2.96 s. 

We interpret these estimates of the rms diffences as representing the average overall errors in 
the A, B, and C grade measurements. (Unmodelled lower mantle structure should be nearly 
identical for data with bounce points within 80 km and at similar azimuths). Under this 
interpretation we can write, for A-quality residuals. 


-o'? 2 y > 

<JArms Z = °Ain + °epi 

Substituting values of o Ams and o epi into (A2, yields o A ni= 0.87 s. Similarly, for B and C quality 
measurements, we find o Bm = 194 s and a Cm = 2.86 s. From (Al). the total unceriainty for A, B, 

and C quality is, respectively, = 1-25 s, Og = ^-14 s ’ s ‘ 

In the weighted regression experiments the A, B, and C quality measurements are weighted 

inversely by their measurement variance. 
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APPENDIX B: ERRORS IN THE ALONG-AXIS 


PROFILES AND CONSTRUCTION OF THE DATA 


COVARIANCE MATRIX 


Errors in Bathymetry, Geoid, and Travel Tin* Profiles 

Uncertainties in the along-axis ptofiies of geoid, bathymeuy, and navel nmes are tmportan 

infection in the inversion. The gridded bathymetiic data , W. M- Oce^Uc Office 
,985) include correcdons for the deviation of water column acoustic velocity fan the assumed 
value of 1500 m S' 1 . The geoid dau, provided in die form of a 0.25* x 0.25* grid [Marsh et at., 
,986], include conections for orb., enors, instrument and atmospheric propagation effects, an 

solid Earth and ocean tides. 

We have averaged the bathymeny and geoid heigh, values within a f x box centered a, 
each SS bounce point. The averaging yields a repmsentative value for a region over app—ly 
one horizontal wavelengti, of die SS wave and acts to smooth on, short-wave, engtb vacations. 

Both bathymetry and geoid are corrected for subsidence arid, seafloor age, using the plate coo mg 
model parsons a^Scla.er, ,977; Porsors a^H^er, ,980,. Error intioduced into depth and 
geoid anomalies by isochron mislocation is difficult to estimate precisely, but for an error in age of 
2 My, depth and geoid errors at 80 My would be about 30 m and 0.2 m, respectively, whde at ^ 
My an error in age of 2 My would have a much larger affect, giving depth and geoid errors o 
m and 0.3 m. The magnitude of this error highlights die imporiance of accurate age determmauon, 

especially at young ages. 

The presence of oceanic sediments is another source of error. In the Adanuc Ocean, die 
sediment thickness increases regularly from less than 100 m along the Mid-Atlantic Ridge toward 
continental margins where it can exceed 1 hm [Ewing e, of, 1973; TucHoUce. 1986], A Mm 
sediment thichness leads to conections to residua, depth and geoid of about 500 m and 0.3 m 
respectively [Cnzenove e, at., ,988; Sheehan and McNun. 1989,. On Atlantic lithosphem of ,00 
My age or less the sediment thichness is less titan 500 m in mo* areas. Hence, neglecting the 
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sediment loading collection should not be ctucial in this region. 

The along-axis profile of SS-S residual is a weighted moving average of 10 adiacent data 
points grouped by latitude, using the weights discussed in Appendix A. The same weights and 
moving average are applied to geoid and bathymetry values at a given SS bounce point, even 
though bathymetry and geoid data are presumed to be of equal quality, in order that these profiles 
will be consistent with the SS-S residuals. The standtud error of the mean values for SS-S residual 
tanges from 0.2 s to 1 .6 s. For bathymetry the range of standard deviations from the mean value is 
from 24 to 370 m, and for geoid. 0.08 to 1.0 m. The largest variances in the bathymetry and geo.d 
data come from the Iceland region (nor* of 60’N), and may be due to the more complicated 

tectonics of this region [White, 1988]. 

Before Fourier transforming, the along-axis profiles must be interpolated to a constant 

spacing We use a simple linear inttnpolation scheme to estimate values at a 0.5- spacing. We 

estimate that the typical error in the interpolated da* is comparable to that in the along-axis moving 
averages, which for bathymeriy is on the order of .25 m. for geoid 0.4 m, and for navel time 1 s. 


Effect of Crustal Thickness Variations 

Our poor fit to topography in the inversion experiments can be at least partially attri 
unmodelled effects such as crustal thickness differences. Variations in oceanic crustal thickness 
about die typical value of 6-7 km [Spudich and Orcutt, 1980; White, 1984; Purdy and Detrick, 

1986) are generally though, to be small a, horizontal scales of 100 km and greater. However, the 
crust beneath the Azores plateau is estimated to be between 8 and 9 km thick [Sear/e, 1976, 
Whitmarsh e.ai, 19821 and that beneath Iceland is a, leas. 8 to 14 km thick [Bjornsscn, 1983). By 
simple isostatic mass balance, the depth anomaly due to excess crustal thickness in the Azores 
region would be about 400 m, and a, Iceland, 200 m to 1.6 km. In genetal, simple variations in 
crustal thickness are insufficient to produce a significant SS-S residual. For crustal and mantle S 
wave velocities of 3.5 and 4.4 km/s, a 2-km variation in crustal thickness would contribute less 
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„» 0 2 s «o an SS-S differential tiavel time cohered for differences in bathymetry. However, a. 
Iceland, where rhe crus, is estinured to be as much as 14 km tirick, rhe additional SS-S tiavel rune 

could be up to 0.8 s. 


Data Covariance Matrix 

The data covariance matrix Rdd IS * c ^ orm 


Rdd = 
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where a„2, On 2 , and o, 2 are the nominal variances of the barhymerry, geord, and travel tune 
data, respectively. We may choose ,o consmrct the data covariance matrix not ro reflect the tine 
variance of the dam but tamer to aliow weighting fc.ween the diffeten. dam sets. In tins way, the 
data covariance matrix can be altered to test me telative contributions of diffemn, dam sets to me 

inversion results. 

In all of our inversions, the covariance matrix is constructed to weight me three dam 
approximately equally. For example, examination of Figure 10 indicates that at 3000- 
wavelength me amplitude of me geoid signal is approximately 4 m, bathymetry 1 km, and tiavel 
time 2 s Thus if a value of 1 m is chosen for o N , then a value of 0.5 s for a, and 0.25 km for o h 
should yield approximately equal weighting of dam sets. The conesponding I/O* values ate men 1 

for geoid, 4 for travel time, and 16 for bathymetry. 
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Figure 1. 


Figure 2. 


Figure 3. 


Figure Captions 


An example of the measurement of SS-S differential travel time for the event 
of December 24, 1985 (10 km focal depth), at GDH (63* epicentral 

distance), (a) “Synthetic” SS pulse generated from S. The S pulse is 

windowed and attenuated to account for the additional time that SS travels in 
the mantle (t* = 3 s), and a rc/2 phase shift is applied, (b) Windowed SS 
wave pulse, (c) Cross-correlation of the trace in (b) with that in (a). The 
differential travel time residual is -5.04 s. 

Distribution of earthquakes (triangles) and seismograph stations (circles) 
used to measure SS-S differential travel times. Stations are from the 
GDSN, NARS, and GEOSCOPE digital arrays. Earthquakes are from the 
Harvaid CMT catalogue (generally mb > 5.0) from the years 1977-1987. 
Lambert equal area projection with pole of projection at 45 N, 40 W. 

(a) Map view of SS-S residuals relative to PREM [Dziewonski and 
Anderson, 1981], corrected for Earth ellipticity and seafloor bathymetry. 
Residuals are plotted at the SS bounce point The size of each symbol 
scales linearly with magnitude of the residual. Lambert equal area 
projection with pole of projection at 40*N, 60* W. Negative residuals 
indicate either early SS or late S. Plate boundaries are from DeMets et al. 

[1990]. 

(b) Same as (a) but including data only for SS bouncepoints on lithosphere 
younger than 100 My. 
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Figure 4. 


Figure 5. 


Figure 6. 


Figure 7. 


Figure 8. 


SS-S travel time residual versus square root of seafloor age for data from 0- 
60*N. Each plotted point represents the weighted mean of 14 adjacent data 
points. Weights are constructed from variances determined as discussed in 
Appendix A. Horizontal and vertical bars are standard errors of the means 
of the travel time residuals and (age)W. Linear regression yields a slope of 
-0.68 ± 0.08 s / (My)W for a 0-100 My age range (solid line) or -0.76 ± 

0.09 s/ (My) 1/2 for a 0-80 My range (dashed line). 

Age-corrected SS-S residual (see text) versus azimuth 0. Each plotted point 
represents the weighted mean of 10 adjacent data points. The solid curve 
shows the best-fitting 40 variation derived from these data. The dashed 
curve shows the preferred model of Kuo et a/. [1987], which corresponds 
to an alignment of the a axis of olivine in the approximate direction N13 W. 

(a) Map view of age-corrected SS-S residuals. 

(b) Same as (a) but including data only with SS bounce points on 
lithosphere younger than 100 My. 

Map view of the distribution of sampling azimuths. Lines indicate the wave 
path azimuth at the SS bounce point. Mercator projection. 

Age-corrected SS-S residual versus latitude along the Mid-Atlantic Ridge 
from 10 to 90*N. The residuals shown are moving averages (such that each 
point is used twice) of 12 adjacent data points. Bounce points on 
lithosphere of age 0-100 My are used. The approximate locations of several 
fracture zones (Fifteen-Twenty, Kane, Atlantis, Oceanographer, and 
Charlie-Gibbs, denoted by abbreviations) and of the Iceland and Azores 
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Figure 9. 


Figure 10. 


Figure 11. 


hotspots are indicated 


Linear correlation, by highest harmonic degree removed from the geoid, of 
observed SS-S residual with geoid height measured at the corresponding SS 
bounce point. Both travel time and geoid residuals are age-corrected. First 
the raw [Marsh et al . , 1986] geoid data are correlated with SS-S residuals 
and a slope and correlation coefficient determined. Then a geoid reference 
field [Lerch et al., 1979] up to degree and order 2 (with taper to degree and 
order 6) is calculated and removed from the geoid data, the slope and 
correlation coefficient with SS-S calculated, and so on for higher harmonic 
degrees (removed from the geoid data, with appropriate tapers (up to (+ 4). 
(a) Linear correlation coefficient between geoid and SS-S residuals vs. 
highest harmonic degree and order removed from the geoid. (b) Slope of 
the correlation between geoid and SS-S data, as a function of highest 
harmonic degree and order removed from the geoid. Extra points at degree 
and order 10 are obtained by using different tapers (no taper, taper to (= 14, 
and taper to (= 15). 

Comparative plots of age-corrected (a) bathymetry, (b) geoid, and (c) SS-S 
residual along the Mid-Atlantic Ridge, 10-65*N. Bathymetry and geoid 
have been high-pass filtered (see text). All of the residuals shown are 
moving averages of 10 adjacent data points. Bounce points from 
lithosphere of age 0- 100 My are used, except that data from the Labrador 

Sea region are omitted. 

Upper mantle portion of the kernels for geoid and topography at two 
wavelengths X = 2rt/k for two viscosity models. The convecting region in 
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Figure 12. 


Figure 13. 


both models is overlain by a high-viscosity layer 40 km thick, with 
viscosity 10 4 that of the underlying mantle. 

(a) High-viscosity lid is underlain by a mantle of uniform viscosity and 
other physical parameters. 

(b) High-viscosity lid is underlain by a zone extending to a depth of 200 km 
having a viscosity equal to 0.01 that of the underlying mantle. 

Upper mantle portion of the kernels for geoid and topography at two 
wavelengths in spherical versus cartesian coordinates. The convecting 
region is overlain by a high-viscosity layer 40 km thick in each of two 
models for viscosity structure. 

(a) Underlying mantle is of uniform viscosity. Cartesian kernels are for 
4000 km wavelength (solid lines) and spherical kernels are for C- 10 
(dashed lines). 

(b) High-viscosity lid is underlain by zone extending to 200 km depth 
having a viscosity equal to 0.01 that of the underlying mantle. Cartesian 
kernels are for 4000 km wavelength and spherical kernels are for (= 10. 

(c) Same as (a) but with Cartesian kernels for 6667 km wavelength and 
spherical kernels for (= 6. 

(d) Same as (b) but with Cartesian kernels for 6667 km wavelength and 
spherical kernels for C= 6. 


Results of combined inversion of geoid, bathymetry, and SS-S travel time 
residuals for upper mantle temperature variations. The viscosity structure is 
taken to consist of a thick high-viscosity lid overlying a constant-viscosity 


halfspace. 

(a) Three solutions for along-axis temperature variations: Dotted line: 
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Figure 14. 


Figure 15. 


Temperature perturbations constrained to be uniform over 0-150 km depth. 
Long-dashed line: Temperature perturbations constrained to be uniform 
over 0-300 km depth. Short-dashed line: Temperature perturbations 
constrained to be uniform over 0-650 km depth. 

(b) Observed (solid line) and predicted along-axis profiles of SS-S travel 
time residual. The “observed” profile is actually a filtered version of the 
observations, containing only the wavelengths used in the inversion (1400 
to 7 100 km). Predicted profiles were calculated from equation 5. Line 
types correspond to those of the temperature models. 

(c) Observed and predicted along-axis geoid profiles. Same treatment as in 

(b). 

(d) Observed and predicted along-axis bathymetry profiles. Same treatment 
as in (b). 

Same as Figure 1 3 except for that the viscosity structure includes a zone 
extending from the base of the lid to a depth of 200 km with a viscosity 
equal to 0.01 that of the underlying mantle. 

Results of combined inversion of geoid, bathymetry, and SS-S travel time 
residuals for variations in upper mantle composition (Mg#). The viscosity 
structure is taken to consist of a thick high-viscosity lid overlying a 
constant-viscosity halfspace. 

(a) Three solutions for along-axis composition variations: Dotted line. 
Composition perturbations constrained to be uniform over 0-150 km depth. 
Long-dashed line: Composition perturbations constrained to be uniform 
over 0-300 km depth. Short-dashed line: Composition perturbations 
constrained to be uniform over 0-650 km depth. 
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Figure 16. 


Figure 17. 


(b) Observed (solid line) and predicted along-axis profiles of SS-S travel 
time residual. 

(c) Observed and predicted along-axis geoid profiles. 

(d) Observed and predicted along-axis bathymetry profiles. 

Same as Figure 15 except for that the viscosity structure includes a zone 
extending from the base of the lid to a depth of 200 km with a viscosity 
equal to 0.01 that of the underlying mantle. 

Results of combined inversion of geoid, bathymetry, and SS-S travel time 
residuals for both upper mantle temperature and composition variations. 

The viscosity structure is taken to consist of a thick high-viscosity lid 
overlying a constant-viscosity halfspace. 

(a) Three solutions for along-axis temperature variations: Dotted line: 
Composition perturbations constrained to be uniform over 0- 150 km depth. 
Long-dashed line: Composition perturbations constrained to be uniform 
over 0-300 km depth. Short-dashed line: Composition perturbations 

constrained to be uniform over 0-650 km depth. 

(b) Three solutions for along-axis composition variations: Dotted line: 
Composition perturbations constrained to be uniform over 0-150 km depth. 
Long-dashed line: Composition perturbations constrained to be uniform 
over 0-300 km depth. Short-dashed line: Composition perturbations 
constrained to be uniform over 0-650 km depth. 

(c) Observed (solid line) and predicted along-axis profiles of SS-S travel 
time residual. 

(d) Observed and predicted along-axis geoid profiles. 

(e) Observed and predicted along-axis bathymetry profiles. 
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Figure 18. 


Figure 19. 


Figure 20. 


Same as Figure 17 except for that the viscosity structure includes zone 
extending from the base of the lid to a depth of 200 km with a viscosity 
equal to 0.01 that of the underlying mande. 

Same as Figure 17 but compositional variations constrained to be from 0-50 
km depth only. 

Same as Figure 18 but compositional variations constrained to be from 0-50 


km depth only. 



TABLE 1. Digital Seismograph Stations Used 


Station Code Network 


ALQ 

DWWSSN 

ANMO 

SRO 

ANTO 

SRO 

BCAO 

SRO 

BER 

DWWSSN 

BOCO 

SRO 

COL 

DWWSSN 

GAC 

CAN 

CDH 

DWWSSN 

GRFO 

SRO 

JAS1 

DWWSSN 

KBS 

DWWSSN 

KEV 

DWWSSN 

KON 

HGLP 

KONO 

ASRO 

UDN 

DWWSSN 

NE04 

NARS 

NB06 

NARS 

NE09 

NARS 

NE10 

NARS 

NE11 

NARS 

NE12 

NARS 


Latitude (°N) Longitude (°E) 


34.942 

-106.458 

34.946 

-106.457 

39.869 

32.794 

4.434 

18.535 

60.387 

5.326 

4.587 

-74.043 

64.900 

-147.793 

45.70 

- 75.47 

69.250 

-53.533 

49.692 

11.222 

37.947 

-120.438 

78.917 

11.924 

69.755 

27.007 

59.649 

9.598 

59.649 

9.598 

46.750 

-121.810 

52.810 

6.670 

50.100 

4.600 

44.850 

0.980 

43.090 

-0.700 

41.480 

-1.730 

40.640 

-4.160 



NE13 


NARS 


38.690 


-4.090 


NE14 

NARS 

NE15 

NARS 

NE16 

NARS 

NE17 

NARS 

RSCP 

RSTN 

RSNT 

RSTN 

RSNY 

RSTN 

RSON 

RSTN 

RSSD 

RSTN 

SCP 

DWWSSN 

SSB 

GEOSCOPE 

TOL 

DWWSSN 

WFM 

GEOSCOPE 


37.190 

-3.600 

50.810 

5.780 

45.763 

3.103 

39.881 

-4.049 

35.600 

-85.569 

62.480 

-114.592 

44.548 

-74.530 

50.859 

-93.702 

44.120 

-104.036 

40.795 

-77.865 

45.280 

4.540 

39.881 

-4.049 

42.610 

-71.490 


ZOBO 


ASRO 


-16.270 


-68.125 



TABLE 2. Adopted Constants 


Variable 

Description 

Value 

a 

volumetric coefficient of thermal 
expansion 

2.5 x 10-5 K' 1 (a) 

PO 

average mantle density 

3300 kg m 

pw 

density of seawater 

1000 kg m° 

6 67 x 10" 11 N m 2 kg' 

r 

gravitational constant 


g 

avs/9T 

surface gravitational acceleration 

9.8 m s' z 

-6.0 x 10' 4 km/s K' 1 

thermal coefficient of 
shear velocity 

dvs/3Mg 

variation of shear velocity with Mg# 

1.8 x lO- 2 km/s/Mg# 

9p/9Mg 

variation of density with Mg# 

-12 kg/m 3 /Mg# (a) (b) 

P 

average SS ray parameter at 70 

0.1375 s/km 


(a) Stacey [1977], Duffy and Anderson [1989] 

(b) Akimoto [1972] 



TABLE 3. Inversion Models and Variance Reduction 


Model: Temperature variations only 
Layer Viscosity 


Variance reduction, % 


thickness 

0-150 km 

structure 

cvm 

m range 

180 K 
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25 

79 

58 

0-150 km 

lvz 

230 K 

57 

27 

79 

66 

0-300 km 

cvm 

60 K 

47 

21 

85 

41 

0-300 km 

lvz 

110 K 

57 

24 

85 

65 

0-650 km 

cvm 

20 K 

41 

14 

91 

25 

0-650 km 

lvz 

33 K 

49 

17 

83 

51 

Model: Compositional variations only 

Variance reduction, % 

Layer 

Viscosity 

structure 


total 

topo 

geoid , 


0-150 km 

cvm 

1.1 

33 

46 

74 

-9 

0-150 km 

lvz 

2.4 

44 

75 

76 

-9 

0-300 km 

cvm 

0.4 

33 

29 

87 

-6 

0-300 km 

lvz 

1.3 

43 

73 

73 

-7 

0-650 km 

cvm 

0.1 

32 

19 

93 

-3 

0-650 km 

lvz 

0.5 

49 

65 

86 

+5 

Model: Thermal and compositional variations in 

same layer 

Variance reduction, % 

Layer 

Viscosity 

structure 

AT 

ranee 

AMg# 
ranee 

total 

tOPO 

geoid 

0-150 km 

cvm 

210 K 

1.5 

75 

44 

78 

0-150 km 

lvz 

235 K 

2.1 

86 

75 

80 


100 

100 



0-300 km 

cvm 

110 K 

0-300 km 

lvz 

125 K 

0-650 km 

cvm 

55 K 

0-650 km 

lvz 

60 K 


0.7 

73 

28 

89 

100 

1.1 

84 

74 

76 

100 

0.4 

71 

18 

94 

100 

0.8 

85 

66 

88 

100 


Model: Thermal inversion in 


layers as noted, compositional variations 0-50 km only 



0-300 km 
0-650 km 
0-650 km 


lvz 

cvm 

lvz 


120 K 
25 K 
35 K 


4.7 

6.0 

4.6 


86 

73 

75 


90 

85 

82 


77 

92 

73 


89 

47 

71 


cvm = constant viscosity mantle 
lvz = mantle with low viscosity zone 
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